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ABSTRACT 


A  new  concept  called  here  the  "trace  function"  was  formu¬ 
lated  and  developed  in  this  work. 

The  notion  is  that  the  phase  differences/time  delays  of 
the  signals  received  from  an  array  of  sensors  which  has  a  given 
defined  geometry  describe  a  certain  "trace  function"  which  is 
constant  for  that  array  at  a  given  look  direction.  If  one  can 
find  a  geometry  for  which  the  trace  function's  basic  shape  is 
not  dependent  on  the  look  direction,  but  only  some  of  its 
parameters  are,  then  one  has  a  powerful  method  to  correlate 
(compare)  the  trace  function  of  the  received  signal  to  a  stored 
replica . 

The  concept  of  the  trace  function  was  formalized  and  its 
application  to  some  typical  array  geometries  was  demonstrated. 
Furthermore,  it  was  shown  that  for  highly  symmetric  arrays  like 
the  circular  and  linear  arrays  the  trace  function  reduces  to  a 
particularly  simple  form. 

This  characteristic  was  used  to  derive  a  relatively  simple 
and  manageable  MMSE  estimator  for  the  bearing  of  the  incoming 
signal.  The  estimator  is  applicable  to  either  narrow  band 
signals  by  use  of  phase  difference  trace  functions  or  to  a 
wide  band  signal  using  the  time  delay  trace  function. 

The  performance  of  the  estimator  was  checked  by  simulation 
and  compared  to  the  CRLB  as  adapted  to  this  application. 

Finally,  a  system  configuration  applying  trace  function 
principles  was  outlined  and  the  major  problem  areas  caused  by 
the  specific  application  were  identified,  reviewed  and  some 
suggestions  to  solutions  were  made.  The  principle  can  be 
applied  to  any  other  situation  in  which  phase  differences 
and/or  time  delays  within  an  array  of  sensors  can  be  measured. 
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I .  INTRODUCTION 


A.  PURPOSE 

The  estimation  of  the  bearing ,  frequency  and  amplitude 
of  acoustic  signals  for  passive  sonar  is  normally  accomplished 
with  an  array  of  sensors  which  has  a  defined  geometrical 
arrangement.  In  most  cases  there  are  several  factors  which 
influence  the  geometry,  of  which  the  most  prominent  are  the 
size,  number  of  elements  and  frequency  coverage,  versus  the 
desired  performance.  As  a  result  of  various  tradeoffs,  tra¬ 
ditionally  three  main  geometrical  configurations  have  been 
applied : 

the  linear  array 

-  the  circular  or  spherical  array 

-  the  conformal  array . 

The  linear  array  is  typically  applied  in  fixed  sites  and  in 
towed  configurations  and  as  such  the  size  (length)  is  not  the 
dominant  limitation.  The  circular  or  spherical  array  and  the 
conformal  array  are  applied  normally  in  moving  platforms  and 
the  size  limitations  (maximum  aperture)  is  a  major  factor  in 
the  design.  The  size  limitation  implies  a  limitation  of  an¬ 
gular  accuracy  and  resolution  at  low  frequencies .  The  extent 
of  the  performance  limitation  can  be  reduced  by  proper  pro¬ 
cessing  of  the  array  outputs. 

To  achieve  the  desired  performance,  the  processing  of 
array  data  is  done  in  two  domains:  space  domain  and  time 


domain  and  their  derived  (by  transform)  versions;  spatial 
and  temporal  frequency. 
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Several  approaches  have  been  taken  to  the  above  pro¬ 
cessing,  from  the  simplest  one  which  is  the  delay  and  sum 
beamformer  [35,  37,  33]  through  the  correlation  processor 
[17,  21]  to  the  optimal  processing  [  l,  2  ,  5  »  13*  20>  24>  39] 
with  adaptive  version  as  the  "ultimate"  processor  [12,  20,  42]. 
All  those  techniques  solve  to  some  extent  the  processing  pro¬ 
blem  of  arrays  whose  size  is  such  that  the  interelt.  -  ;nt 
spacing  is  on  the  order  of  X/2  for  the  lowest  frequencies  of 
interest,  and  have  enough  elements  to  sample  adequately  the 
space/time  field. 

Attempts  to  overcome  the  spatial  resolution  problem  have 
been  made  by  applying  the  "superdirective  array"  concept. 
Although  this  technique  gives  infinitely  good  spatial  resolu¬ 
tion  for  noise  free  environment,  it  has  not  been  found  to  be 
practical  because  of  poor  aperture  efficiency  and  sensitivity 
of  the  solution  to  tolerances  in  implementation  [11,  35,  37]. 
Recently  some  constrained  sensitivity  solutions  to  the  super¬ 
directive  arrays  were  presented  which  give  a  less  sensitive 
and  more  efficient  array  [io ,  24], 

Some  of  the  above  processors  are  implemented  as  "beam- 
formers"  while  the  others  are  considered  as  estimators  of  the 
incoming  signal  direction.  As  a  general  conclusion  it  can  be 
said  that  an  adequate  solution  to  the  problem  of  small  arrays 
at  low  frequencies  has  not  yet  been  fully  derived.  The  aim 
of  this  dissertation  was  to  develop  such  a  solution. 
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B.  MOTIVATION 


The  practical  situation  which  motivated  the  initiation 
of  this  work  was  the  small  circular  and  cylindrical  arrays 
found  in  a  great  number  of  submarines  and  surface  ships. 

These  arrays  are  ccnstrained  in  size  to  few  meters  diameter 
which  in  turn  limits  their  ability  to  perform  spatial  pro¬ 
cessing  of  low  frequency  signals.  The  size  of  the  whole 
array  is  small  compared  to  the  wavelength  at  frequencies  of 
100  Hz  or  less. 

However,  those  low  frequencies  are  the  very  ones  which 
offer  longer  detection  ranges  and  better  classification  means. 
The  spatial  processing  problem  can  be  solved  by  adding  a  long 
line  array  (towed  in  general)  which  will  give  the  needed 
aperture  for  the  low  frequency  spatial  resolution.  For  small 
platforms  this  approach  may  not  be  practical  and  a  different 
way  has  to  be  found  by  signal  processing  of  the  small  circu¬ 
lar  array . 

C.  OUTLINE  OF  THE  DISSERTATION 

A  novel  idea  which  is  believed  to  be  original  is  used  as 
the  basis  for  this  dissertation.  This  idea  is  basically  a 
new  way  to  look  at  the  available  information  from  an  array 
of  sensors  called  the  "trace  function"  by  the  author. 

It  is  based  on  the  fact  that  the  directional  information 
of  an  incoming  plane  wave  is  mainly  contained  in  the  received 
signal  phases,  or  to  be  more  precise,  in  the  relative  phase 
differences  between  received  signals  at  the  various  elements 
of  the  array. 
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Those  phase  differences  are  representative  of  the  time 
delay  of  the  plane  wave  between  the  elements  of  the  array. 

The  notion  is  that  the  phase  differences/time  delays  of  the 
signals  received  from  an  array  of  sensors  which  has  a  given 
defined  geometry  describes  a  certain  "trace  function"  which 
is  constant  for  that  array  at  a  given  look  direction.  If 
one  can  find  a  geometry  for  which  the  trace  function's  basic 
shape  is  not  dependent  on  the  look  direction,  but  only  some 
of  its  parameters  are,  then  one  has  a  powerful  method  to 
correlate  (compare)  the  trace  function  of  the  received  signal 
to  a  stored  replica. 

Fortunately  one  such  geometry  is  provided  by  the  circular 
array,  which  was  the  motivation  behind  this  work  and  whose 
phase  difference  trace  function  is  a  simple  two  dimensional 
sine-cosine  "wave"  of  exactly  one  period  of  each  dimension. 

The  apparent  "phase"  of  this  trace  function  is  indicative 
of  the  bearing  of  the  incoming  plane  wave  signal.  The  trace 
function  idea  can  be  used  with  array  geometries  other  than 
the  circular  one,  however,  it  will  be  less  efficient  as  the 
symmetry  is  broken.  For  example,  it  will  be  shown  in  later 
chapters  how  it  applies  to  linear  arrays. 

General  Definition  -  A  trace  function  is  a  discrete  plot 
of  phase  difference/time  delay  versus  element  indices  for  all 
element  pairs. 

Fig.  1-1  shows  such  a  trace  function  for  a  random  array 
while  Figs.  1-2  and  (-3  show  one  each  for  a  linear  and 
circular  array.  ¥ 
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DIFFERENCE  TRACE  FUNCTION  FOR  fl  LINEAR  ARRAY 


PHASE  DIFFERENCE  TRACE  FUNCTION  FOR  R  CIRCULAR  ARRAY 


Figure 


Two  main  features  can  be  readily  observed  for  any  trace 
function: 

-  the  elements  on  the  diagonal  are  always  zero  (At^=0 
and  A<f>^=0). 

-  the  trace  function  is  always  antisymmetric  across  the 

diagonal  ( Ar .  .  =  -Ax.  .  and  Ad> .  .  =-A<|>.  . ) . 
s  i]  31  T3i 

For  the  noise  free  conditions  the  trace  function  can  be 
derived  easily  for  a  given  geometry.  Trace  functions  are 
defined  separately  for  each  look  direction  and  for  each  fre¬ 
quency  resolution  cell  (or  for  wide  band  time  delay). 

In  Chapter  II  of  this  dissertation  the  subject  of  defini¬ 
tion  and  construction  of  the  trace  function  will  be  treated 
in  great  detail. 

For  the  received  signal  in  a  noisy  environment  the  above 
trace  functions  must  be  estimated  for  each  narrow  band  signal 
present  and  for  the  wide  band  data. 

Once  the  trace  function  of  the  received  signals  is 
estimated  it  can  be  "matched"  to  a  stored  replica  of  the 
noise-free  trace  function  of  the  array  (either  in  discrete 
or  analytic  form) . 

The  task  of  estimating  the  trace  function  from  the  re¬ 
ceived  data’  is  not  trivial  mainly  because  it  has  to  be  done 
in  the  very  difficult  conditions  of  small  array  size  and 
highly  correlated  noise  between  sensors. 

All  optimal  processors  use  some  scheme  to  "whiten"  the 
noise  in  spatial  domain  and/or  temporal  domain  before  it  is 
further  processed.  This  approach  has  to  be  considered  in 


T 


this  case  too.  However,  the  trace  function  concept  poses  a 
heavy  demand  on  that  type  of  preprocessing  because  the 
preservation  of  relative  phase  within  the  array  is  required. 
Present  prewhitening  schemes  are  not  readily  applicable  and 
ways  have  to  be  found  to  perform  the  above  task.  Suggestions 
about  relative  phase  preserving  spatial  decorrelation  methods 
are  made  in  Chapter  V  (System  Considerations)  which  use  a 
compensation  scheme  to  preserve  phase,  and  are  based  on  the 
knowledge  (a  priori  or  estimated)  of  the  ambient  noise  spec¬ 
trum  and  covariance  functions. 

The  estimation  of  phase  differences  can  be  done  by  es¬ 
tablished  cross  spectral  estimation  methods,  especially  the 
"short  modified  periodogram"  method  [40,  25,  26]  which  was 
adopted  in  this  work  and  is  reviewed  in  Chapter  V  for  this 
application. 

For  the  time  delay  estimation  several  approaches  are 
possible  and  this  subject  is  treated  extensively  in  the 
current  literature  [6,  9,  14,  15]. 

Of  the  above  methods,  the  one  which  was  the  most  appli¬ 
cable  is  reviewed  in  Chapter  V  and  some  adaptations  to  this 
work  are  suggested . 

The  match  of  the  estimated  and  stored  trace  function  was 
done  by  a  minimum  mean  square  error  (MMSE)  estimator  which  was 
fully  derived  in  this  work.  The  derivation  of  the  estimation 
is  given  in  Chapter  III  and  its  performance  is  detailed  in 
Chapter  IV.  The  performance  of  the  estimator  was  checked  by 
computer  simulations  and  the  results  are  presented  and  compared 
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to  a  classical  benchmark,  the  Cramer-Rao  Lower  Bound  (CRLB). 

A  block  diagram  for  a  possible  implementation  scheme  is 
given  in  Chapter  V. 

Chapter  VI  presents  a  concise  summary  of  this  work.  In 
addition  it  is  suggested  that  the  trace  function  concept  and 
the  estimators  derived  in  this  work  are  not  limited  in  their 
application  to  sonar.  They  can  be  applied  to  related  areas 
like  seismic,  and  ocean  wave  analysis.  The  chapter  concludes 
with  some  proposed  topics  for  future  work  which  were  opened 
up  by  the  introduction  of  the  trace  function  principle. 

D.  RELATED  WORK 

As  the  concept  is  new,  very  little  related  work  can  be 
mentioned  as  applying  strictly  to  the  basic  idea  of  the  trace 
function.  Only  two  works  [3,  221  can  be  defined  as  such. 

The  rest  of  the  references  used  in  this  dissertation  could  be 
more  correctly  described  as  supporting  works. 

Two  previous  papers  published  by  the  author  (with  G.  L. 
Sackman)  [30,  31]  on  the  subject  of  the  dissertation  are  also 
cited. 

The  two  paper's  which  are  related  to  the  subject  [3,  2?, 
op.cit.]  were  brought  to  the  attention  of  the  author  recently 
after  the  trace  function  concept  had  already  been  developed. 
The  relationship  is  restricted  to  the  general  concept  of 
"matching"  some  function  of  the  cross  spectral  matrix  to  an 
incoming  plane  wave.  Those  two  papers  were  analyzed  in  order 
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to  contrast  with  and  emphasize  the  present  work. 

The  first  paper  by  Munk  et  al  [22]  is  the  basic  one  in 
which  the  authors  address  the  problem  of  estimation  of  the 
direction  of  ocean  waves  arising  from  distant  storms  using  a 
relatively  small  array  of  three  sensors.  The  problem  is 
somewhat  similar  to  the  present  one,  however  it  is  in  a 
totally  different  discipline  and  the  methods  used  for  solu¬ 
tion  are  also  significantly  different.  The  solution  which 
was  shown  in  that  paper  is  a  MSEE  fitting  of  a  full, 
measured  cross  spectral  matrix  (complex  in  general)  of  an 
array,  to  an  assumed  form  of  such  a  matrix  for  a  plane  wave. 
Their  method  used  an  "educated"  search  technique  to  minimize 
the  error  in  the  fitting  process. 

Very  good  results  were  achieved,  however  the  computa¬ 
tional  complexity  was  tremendous. 

The  second  paper  by  Bennett [  3  ]  is  a  tutorial  report 
published  almost  ten  years  after  the  first  one.  In  this 
report  the  author  formalizes  and  reviews  the  methods  of  the 
paper  by  Munk  et  al  adding  some  alternatives  and  demonstrating 
the  theory  by  the  results  of  processing  real  data.  However, 
no  departure  was  made  from  the  basic  approach  of  Munk  et  al. 

It  is  worthwhile  to  mention  that  Bennett  concluded  that 
this  method  of  least  square  fitting  was  the  best  available 
at  that  time  for  bearing  estimation.  In  contrast  to  the 
above  the  present  work  uses  only  certain  relevant  parts  of 
the  cross  spectral  matrix  -  the  phase  differences  -  as  a 
measure  to  be  fitted  and  introduces  the  notion  of  the  trace 
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function  as  a  replica  which  is  only  dependent  on  array  geo¬ 
metry.  This  simplification  leads  to  an  explicit  MMSE 
estimator  equation  for  some  well  defined  array  geometries. 

The  option  of  "educated”  search  remained  open,  but  for  a 
much  simpler  set  of  data.  All  this  is  done  without  degrading 
performance;  in  fact  it  is  demonstrated  by  simulation  that 
the  present  estimator  is  asymptotically  efficient  (achieves 
C.  R.  Lower  Bound). 

Another  fact  which  should  be  mentioned  is  that  the  deri¬ 
vation  of  the  trace  function  concept  was  initiated  from  the 
end  product,  the  bearing,  and  how  it  effects  the  various 
measurable  parameters  in  the  array  and  not  from  the  mathemati 
cal  structure  of  the  cross  spectral  matrix  as  in  the  previous 
works . 

Nevertheless,  the  merit  of  these  related  works  remains 
valid  and  both  the  previous  and  the  author’s  methods  comple¬ 
ment  each  other,  opening  perhaps  a  whole  new  way  of  pro¬ 
cessing  array  data.  The  supporting  works  will  not  be  con¬ 
sidered  in  this  section  mainly  because  they  cover  a  wide 
variety  of  subjects  and  do  not  affect  the  heart  of  the  pro¬ 
blem.  They  are  mentioned  and  referenced  elsewhere  in  this 
work,  each  along  with  its  applicable  subject. 

E.  CONTRIBUTIONS 

The  principal  contributions  of  this  dissertation  can  be 
summarized  as  follows: 

1.  The  first  introduction  and  derivations  of  the  phase- 
difference/time-delay  trace  function  concept  for  arrays  of 
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sensors,  its  applicability  and  limitations. 

2.  Derivation  of  a  MMSE  bearing  estimator  which 
implements  the  trace  function  concept  for  circular  and 
linear  arrays,  and  its  performance. 

3.  Definition  of  a  system  which  applies  the  new  concepts 
along  with  suggestions  for  solution  of  some  specific  problems 
which  arose  from  the  special  demands  of  phase  preserving 
preprocessing  for  the  trace  function  concept. 


II.  PHASE  DIFFERENCE/TIME  DELAY  TRACE 

TOlWTiT "ARRAYS"  "-"THE  gMggPT" 

A.  INTRODUCTION 

The  basic  concept  of  the  trace  function  of  an  array  is 
the  basis  for  the  bearing  estimation  method  developed  in  this 
dissertation. 

In  order  to  give  the  reader  a  good  understanding  of  the 
subject,  the  concept  of  the  trace  function  is  defined  and 
developed  for  a  general  three  dimensional  array  first  and 
then  for  some  typical  two  dimensional  geometries  of  arrays. 
Finally  the  advantages  of  this  concept  will  be  pointed  out 
for  symmetrical  arrays,  especially  the  circular  array. 

Shown  in  Fig.  II-l  is  an  arbitrary  three-dimensional 
array  of  N  receiving  elements.  The  location  ot  the  i-th  ele¬ 
ment  with  respect  to  a  fixed  (x,y,z)  cartesian  coordinate 
system  is  defined  by  the  vector  r^  specifying  both  distance 
and  direction  from  the  center  of  the  coordinate  system. 

A  uniform  plane  wave  propagates  with  speed  c  in  a 
direction  a  with  respect  to  the  coordinate  system.  The 
vector  a  has  a  unit  length. 

The  plane  wave  is  assumed  to  contain  a  wide  band  random 
signal,  S  ^  and  a  set  of  M  line  (sinusoidal)  components  with 
random  initial  amplitude  and  phase,  $nt).  Then,  if  we  assume 
that  the  same  plane  wave  arrives  at  each  element,  we  can 
specify  vector  signal  to  be: 


21 


The  above  definition  of  a  uniform  plane  wave  target  signal 
shows  that  the  signal  component  at  each  element  is  identical 
except  for  a  delay.  This  is  under  the  assumption  that  atten¬ 
uation  is  negligible  and  no  change  in  the  speed  of  sound  occurs 
within  the  array  dimensions 

From  a  total  system  point  of  view  we  have  to  include  also 
the  additive  noise  component  at  each  element;  however,  for 
the  conceptual  demonstration  of  the  trace  function  the  noise 
free  condition  is  adequate.  Under  the  above  stated  conditions 
it  is  suggested  that  these  delays  (t^)  are  the  parameters 
which  contain  the  directional  information  of  the  plane  wave. 

It  is  convenient  to  consider  the  wide  band  portion  of  the 
signal  and  the  line  components  separately,  although  it  will 
be  shown  that  the  trace  function  concept  is  applicable  to  both 
in  a  similar  manner. 

B.  NARROW  BAND  COMPONENT 

As  stated  before ,  the  narrow  band  component  S  ^  is  com¬ 
posed  of  a  set  of  sinusoidal  lines: 

M 

Snb^ji)  =  Z]  cos  [w^(t-r^)  +  (II-3) 

l-l 

where  A^  -  line  amplitude 
<(>£  -  line  phase 

The  phase  difference  trace  function  is  constructed  by 
taking  all  possible  pairs  of  elements  and  observing  that  the 
phase  difference  between  any  pair  of  elements'  signals  Snb(i,H) 
and  S  j  ,  H )  is: 
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<11-4  > 


f 


A$  .  .  .  =  W.  (t  .-T  .  ) 
yi3,£  l  1  ] 


Using  <11-2 )  and  (II-4): 


<rrV 


( 1 1  —  5  ) 


where  "i"  and  "j"  are  equivalent  indices  of  the  array  ele¬ 
ments  and  denote  the  pairs  which  are  considered  for  phase 
difference  or  time  delays. 

Definition:  A  phase  difference  trace  function  is  defined 
as  the  plot  of  A<j>  (phase  difference  as  a  discrete  function  of 
"i"  and  "j"  (element  numbers  in  the  array)  for  "2,"  (frequency) 
constant . 

An  example  of  such  a  trace  function  is  shown  in  Figure 
II-2  for  a  random  array. 

Observations : 

a.  Each  A$  is  a  constant  for  a  given  look  direction 
a  so  that  for  a  fixed  geometry  (fixed  r . )  and  frequency  u). 
the  trace  function  in  invariant  in  time. 

b.  For  various  frequencies  u^=l,2..M  we  get  a  family 
of  such  invariant  trace  functions  which  are  linearly  related. 

c.  For  various  look  directions  we  get  a  set  of  such 
families,  one  for  each  look  direction  a. 

From  the  above  observations  we  can  conclude  that  if  we 
know  such  a  set  of  families  of  trace  functions  which  are 
invariant  for  a  given  array  geometry  we  can  "match"  a  re¬ 
ceived  and  estimated  trace  function  to  the  precomputed  members 
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of  the  sets/f ami lies  and  find  the  look  direction  for  each 
frequency  resolution  element. 

The  task  of  performing  the  above  match  can  be  quite  for¬ 
midable  for  a  practical  situation,  unless  some  regularity  and 
symmetry  can  be  found  within  the  sets.  It  may  even  be  pro¬ 
hibitive  to  be  processed  in  real  time  systems. 

Fortunately,  for  some  practical  and  symmetrical  geometries 
of  arrays,  this  task  of  matching  is  facilitated. 

It  is  easy  to  see  that  this  concept  applies  to  multiple 
targets  as  well  if  we  assume  that  we  can  always  find  a  fine 
enough  frequency  resolution  such  that  two  different  targets 
will  have  their  line  components  in  different  frequency  anal¬ 
ysis  resolution  cells.  (This  condition  is  easily  achieved 
in  practice.)  If  they  are  coming  from  different  look  direc¬ 
tions  they  will  be  members  of  different  families  as  defined 
above  and  they  will  be  recognized  as  such. 


C.  WIDE  BAND  COMPONENT 

For  a  wide  band  random  process  the  phase  is  not  defined, 
hence  we  will  use  the  time  delay  as  a  measure  for  the  trace 
function.  As  mentioned  before  the  wide  band  component  of  the 
signal  is  defined  as  Swt)(t-T^). 

In  analogy  with  the  line  components 


a 


4tij  *  y-y  = 


y-y 


ai-6) 


Definition:  A  time  delay  trace  function  is  defined  as  the 
plot  of  (time  delay)  as  a  discrete  function  of  "i"  and 

"j"  (element  numbers  in  the  array). 
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Figure  II-3  shows  such  a  time  delay  trace  function  for 


a  random  array. 

The  observations  for  the  narrow  band  components  are 
readily  applicable  to  the  wide  band  components  as  well  with 
one  exception:  the  frequency  domain  collapses  into  only  one 
trace  function,  such  that  it  results  in  a  single  set  of  trace 
functions,  one  for  each  look  direction  ct. 

The  same  conclusions  as  for  the  narrow  band  components 
apply,  however  the  multiple  target  distinction  based  on 
frequency  resolution  cannot  be  done  any  more.  Some  sugges¬ 
tions  to  overcome  this  problem  will  be  made  in  Chapter  V. 

D.  TRACE  FUNCTION  CALCULATION 

1.  Analytic  Trace  Functions 

The  trace  functions  of  arrays  which  have  a  geometry 
such  that  they  can  be  expressed  in  a  closed  form  equation  are 
defined  here  as  "analytic".  The  calculation  of  these  trace 
functions  is  simple  and  is  based  on  geometry  only.  In  this 
work  only  two  dimensional  arrays  were  treated,  however  the 
extension  to  three  dimensions  is  straightforward. 

The  calculations  are  derived  from  equations  (II-5) 
and  (II-6)  for  phase  difference  and  time  delay  trace  functions 
respectively. 

These  derivations  are  shown  in  Appendix  A  for  circu¬ 
lar  and  linear  arrays . 

a.  Circular  Array 

The  configuration  of  this  array  is  shown  in  Figure 
II-4 .  From  Appendix  A: 
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-  Time  delay  trace  function 


At . .  =  -2—  sin 
13  c 


[ff  <i+3>-»]sin[j  U-j)] 


( II-7 ) 


-  Phase  difference  trace  function 


iTi 3“» 


2ttD 


sin[l(i+j)  "  9] 


x  sm 


in[j 


(II-8) 


where  R,D 
N 
0 


radius  and  diameter  respectively 
number  of  elements  in  the  array 
incoming  signal  bearing  (look  direction) 


Figure  II-4 

Circular  Array  Configuration 


In  Figure  II-5  through  11-13  are  shown  the  three 
dimensional  plots  of  the  phase  difference  trace  functions  - 
0  for  a  circular  array  of  32  elements  (Equation  II-8) 
at  various  look  directions  (bearings).  The  effect  of  fre¬ 
quency  change  can  be  seen  by  comparing  Figures  II- 5  and  11-14. 
This  effect  is  only  on  the  "amplitude"  of  the  trace  function 
but  not  on  its  basic  form  (similar  effects  are  caused  by 
changes  in  diameter-D  or  sound  velocity-C).  The  computer 
program  for  ploting  these  trace  functions  on  a  Hewlett-Packard 
9845T  computer  is  presented  in  Appendix  B-l. 

Observations : 

(1)  The  trace  functions  of  circular  arrays 
are  always  a  sampled  two  dimensional  combined  cosine  and  sine 
surface  of  exactly  one  period  (360°)  in  each  dimension  ("i" 
and  " j " ) . 

(2)  The  "phase"  of  the  above  two  dimensional 
wave  is  identically  the  direction  of  arrival  of  the  incoming 
signal . 

(3)  The  estimation  of  the  "phase"  can  be 
done  relatively  easily  for  the  noisy  trace  function  by  mini¬ 
mum  mean  square  error  surface  fitting. 

b .  Linear  Array 

The  configuration  of  this  array  is  shown  in 
Figure  11-15. 

From  Appendix  A: 

Time  delay  trace  function 

At..  =  ^(i-j)sinQ  ( II-9 ) 

ij  c 


30 


difference 


DIFFERENCE  TRRCE  FUNCTION  FOR  ft  CIRCULAR  RRRFT 

r  i  nn  —  30  DEG.  Fr@quency=  250  Hz  Di»rrieter=  5 


Figure  II- 


PHASE  DIFFERENCE  TRACE  FUNCTION  FOR  R  CIRCULAR  ARRAY 

Besrino=  45  DEG.  Frequency-  250  Hz  Di  arnate  r  =  5  Meter 


DIFFERENCE  TRACE  FUNCTION  FOR  A  CIRCULAR  ARRAY 


PHASE  DIFFERENCE  TRACE  FUNCT 


DIFFERENCE  TRACE  FUNCTION  FOR  A  CIRCULAR  ARRAY 


Figure  11-10 


PHASE  DIFFERENCE  TRACE  FUNCTION  FOR  A  CIRCULAR  ARRfi' 


Figure  11-11 


FERENC  E  TRACE  FUNCTION  FOR  H  CIRCULAR  ARPH 


Figure  11-12 


PHASE  DIFFERED 


Figure  11-13 


PHASE  DIFFERENCE  TRACE  FUNCTION  FOP  A  CIRCULAR  ARRAY 

Boarir-ig=  0  DEG.  F  rsqusn  cy=  100  Hz  D  i  sriiete  r  =  5  Meter 


i  [  E  1  erne  n  t  #3 
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t-o 

Figure  11-15 

Linear  Array  Configuration 

In  Figures  11-16  through  11-20  are  shown  the  three  dimensional 
plots  of  the  phase  difference  trace  functions  0  for  a 

linear  array  of  32  elements  (Equation  11-10)  of  the  same 
aperture  as  the  circular  array.  The  computer  program  is 
presented  in  Appendix  B-2. 

Observations : 

■  (1)  The  trace  functions  of  a  linear  array  are 

always  a  sampled  two  dimensional  plane  surface. 

(2)  The  "slope"  of  this  plane  is  indicative 
of  the  direction  of  arrival  of  the  incoming  signal. 
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PHASE  DIFFERENCE  TRACE  FUNCTION  FOR  A  LINEAR  ARRAY 

Bearing®  30  DEG.  Frequency®  250  Hz  Length*  5  Meter 


eirie 
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PHASE  DIFFERENCE  TRACE  FUNCTION  FOR  A  LINEAR  ARRAY 

Bearinq*  270  DEG.  Frequency*  250  Hz  Length*  5  Meter 
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(3)  The  estimation  of  this  slope  can  be  done 
easily  for  the  noisy  trace  function  by  minimum  mean  square 
error  surface  fitting. 

2 .  General  Trace  Functions 

For  those  arrays  which  do  not  have  analytic  trace 
functions  the  approach  is  slightly  different.  In  this  case 
in  order  to  use  equations  (II-5)  and  (II-6)  the  value  of  the 
vectors  r^  is  needed  either  in  cartesian.  (x,y)  or  polar 
(r,9)  coordinates.  Then  each  (Eq.  II-2)  is  calculated 
separately  as  referenced  to  the  origin  of  the  coordinate 
system,  and  the  differences  are  calculated  and  plotted.  In 
this  section  two  geometrical  configurations  were  chosen  for 
demonstration: 

a.  A  "Conformal  array"  which  is  composed  of  a  half 
circle  and  two  linear  parts.  The  configuration  is  shown  in 
Figure  11-21.  For  this  array,  although  it  is  built  of  well 
defined  geometrical  forms  its  trace  function  cannot  be 
expressed  in  closed  form.  Figures  11-22  through  11-28  present 
the  trace  function  for  various  look  directions  (bearings). 

It  can  be  observed  that  it  still  preserves  symmetry  around 
the  zero  bearing  and  it  is  built  up  of  pieces  of  the  trace 
functions  of  circular  and  linear  arrays. 

b.  A  "Random  array"  is  illustrated  which  is  composed 
of  19  elements  randomly  distributed  within  a  circle  of  about 

5  meters.  The  configuration  is  shown  in  Figure  11-29.  The 
trace  functions  for  various  look  directions  (bearings)  are 
shown  in  Figures  11-30  through  11-34.  These  trace  functions 


47 


are  totally  random  and  only  two  obvious  types  of  symmetries 
are  preserved: 

-  A*..  =  -A*j;L 

-  the  whole  trace  function  is  antisymmetric  for 
bearings  of  180  degrees  apart  (this  second 
characteristic  is  a  result  of  the  first  one) . 

The  computer  program  which  generates  both  the  "conformal" 
and  "random"  trace  can  be  found  in  Appendix  B-3. 
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PHASE  DIFFERENCE  TRACE 


PHASE  DIFFERENCE  TRACE  FUNCTION  FOR  A  CONFORMAL  ARPA 


PHASE  DIFFERENCE  TRACE  FUNCTION  FOR  A  CONFORMAL  ARRAY 

Bearirici=  180  DEG.  Frequency*  200  Hz 


Figure  11-27 


Figure  11-29 

Random  Array  Configuration 
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PHASE  DIFFERENCE  TRACE  FUNCTION  F< 


Figure  11-30 


PHASE  DIFFERENCE  TRACE  FUNCTION  FOR  A  RANDOM  ARRAY 


Figure  11-31 


PHASE  DIFFERENCE  TRACE  FUNCTION  FOR  A  RANDOM  ARRAY 
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PHASE  DIFFERENCE  TRACE  FUNCTION  FOR  A  RANDOM  ARRAY 
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Common  characteristics  of  the  non-analytic  trace 
functions  are: 

They  do  not  have  a  distinct  single  parameter 
which  is  indicative  of  the  bearing  of  the  incoming  signal; 
however,  they  still  have  a  distinct  form  for  each  bearing. 

As  a  result  of  (a)  the  minimum  mean  square  error 
surface  fitting  still  can  be  used;  however,  each  trace  func¬ 
tion  must  be  fitted  on  a  point-by-point  basis.  Some 
"educated"  search  procedures  may  be  used. 

NOTE:  1.  All  the  trace  functions  plotted  in 

this  part  were  phase  difference  trace 
functions .  The  time  delay  trace 
functions  were  not  given  because 
they  are  only  scaled  versions  with 
the  same  functional  form. 

2.  Each  intersection  in  the  trace  func¬ 
tion  plots  corresponds  to  a  data  point. 

E.  DEGENERATE  TRACE  FUNCTIONS 

In  some  applications  the  task  of  generating  and  processing 
the  full  trace  function  may  not  be  necessary  or  possible.  For 
those  situations  a  degenerate  version  of  the  trace  functions 
can  be  applied.  The  degenerate  .  versions  are  derived  by  taking 
only  some  distinct  pairs  from  all  the  possible  pairs  of  ele¬ 
ments  . 
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For  example  for  a  circular  array  a  logical  choice  of 
pairs  will  be  on  opposite  sides  of  diameters  such  that: 


j  =  (i  +  |)  mod  N  (11-11) 

Then,  by  substitution  in  (II-7) 

At..  =  ^  cos  (2gi  -  0)  (11-12) 

11  C  N 

and  in  (II- 8) 

"  4Tij “i  *  T5  “s  -  6)  (II-13> 

which  is  a  simple  cosine  wave  whose  phase  is  identically  the 
bearing  of  the  incoming  signal.  It  can  be  shown  that  the 
performance  of  the  degenerate  trace  function  is  not  as  good 
as  that  of  the  full  trace  function  but  may  be  adequate  for 
some  applications . 

F .  CHAPTER  SUMMARY 

In  this  chapter  the  concept  of  the  trace  function  was 
developed  and  demonstrated.  From  the  few  examples  shown  it 
can  be  concluded  that  the  concept  is  well  adapted  to  be  used 
in  a  generalized  spatial  replica  correlation  scheme.  In 
addition,  for  highly  symmetrical  arrays  the  amount  of  pro¬ 
cessing  requires  is  significantly  reduced. 

As  such  the  circular  array  is  a  natural  one  to  be  treated 
this  way.  In  subsequent  chapters  the  bearing  estimator  for 
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a  circular  array  will  be  developed  based  on  trace  function 
concepts.  In  addition  to  the  circular  array  the  linear  array 
will  be  also  treated  mainly  for  comparison  purposes. 


III.  MINIMUM  MEAN  SQUARE  ERROR  SURFACE  FITTING  AS 

APPLIED  TO  BEARING  ESTIMATION  BY  TRACE  FUNCTIONS 

A.  INTRODUCTION 

In  the  previous  chapter  it  was  shown  that  the  trace 
function  can  be  represented  by  a  three  dimensional  surface 
(function)  whose  shape  is  dependent  on  the  bearing  0.  For  a 
given  array  the  trace  functions  are  known  either  analytically 
or  in  the  form  of  a  data  matrix  at  any  bearing  0.  For  the 
received  data  from  an  array  of  sensors  the  trace  functions 
are  estimated  via  a  cross  spectral  estimation  procedure,  and 
are  always  in  a  data  matrix  form.  In  order  to  estimate  the 
bearing  by  means  of  the  trace  function  one  has  to  "match" 
the  estimated  trace  function  to  the  known  (stored  or  analytic) 
trace  function.  There  are  several  ways  to  perform  the  above 
"match",  some  of  them  require  the  knowledge  of  the  noise  sta¬ 
tistics  and  others  do  not  require  it.  As  the  noise  statistics, 
especially  the  cross  noise  statistics,  may  be  inaccessible, 
the  methods  which  will  be  followed  in  this  analysis  will  not 
require  this  knowledge . 

The  best  known  method  of  this  class  is  the  minimum  mean 
square  error  (MMSE).  In  this  method  the  sum  of  squared  errors 
between  the  estimated  (measured)  value  of  the  trace  function 
and  the  "known"  value  of  it  at  all  coordinates  "i,j"  will  be 
minimized  with  respect  to  the  parameter  "0"  (the  bearing). 

The  value  of  0  which  will  give  the  minimum  is  the  expected  one. 

In  order  to  give  the  reader  a  better  intuitive  feeling 
about  the  estimation  problem  Figures  III-l  through  IIT-10 
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NOISY  PHRSE  D IFF .  TRRCE  FUNCTION  FOR  fl  CIRCULAR  RPRR 


i  .  >■»  ,'wtt.  Hum*-*  Ml> 


NOISY  PHASE  DIFF.  TRACE  FUNCTION  FOR  A  CIRCULAR  ARRAY 


Figure  III- 


TRRCE 


1  CE  lenient  #3 


NOISY  PHRSE  DIFF.  TRACE  FUNCTION  FOR  fi_ LINEAR  ARRAY 
Bearing3  30  DEG.  L-'^  —  .83  Var=  .001  Rad**2 


-III 


NOISY  PHASE  DIFF.  TRACE  FUNCTION  FOR 

Bearinq=  30  BEIG.  =  .83  Var» 


Figure  III- 


figure 


NOISY  PHASE  DIFF.  TRHCE  FUNCTION  FOR  fl  LINEAR  ARRAN 


Figure  III- 


:SY  PHASE  DIFF.  TRACE  FUNCTION  FOR  A  LINEAR  ARRAY 


Figure  III-10 


show  noisy  trace  functions  one  would  get  after  the  cross 

spectral  estimation  process.  The  above  figures  show  trace 

functions  for  circular  and  linear  arrays  for  various  noise 

2 

variance  levels  from  .001  Rad  (St. dev.  1.86deg.)  to  1.152 
2 

Rad  (St. dev.  61.5deg.).  It  can  be  seen  that  even  for  the 
highest  noise  level  shown  one  can  still  visually  perceive 
the  form  of  the  trace  function  while  for  a  single  data  point 
this  noise  may  be  overwhelming.  The  cause  of  this  perception 
is  the  visual  integration  over  the  surface  of  the  trace 
function.  The  very  same  characteristic  is  used  by  an  esti¬ 
mation  processor. 

B.  DERIVATION  OF  THE  ESTIMATOR 

As  it  was  shown  the  phase  difference/time  delay  trace 
function  C A<|> ^ ^  Q  and  At.^  q)  is  dependent  on  indices  "i,j" 
and  parameter  8 . 

NOTE:  A<() .  .  *  will  be  used  throughout  the  derivation  but 

i]  ,  8 

the  analysis  applies  as  well  to  At..  _. 

1  j  ,  8 

The  estimated  trace  function  will  be: 

Aii-  =  fTr(i,j,e>  ♦  v(i,j)  (III-1) 

where  f^r  is  the  trace  function  and  v(i,j)  is  a  random  esti¬ 
mation  error  which  will  be  approximately  gauss ian  for  long 

enough  averaging  time  in  the  estimation  process.  As  it  was 
* 

shown  in  [26]  the  error  variance  is  not  dependent  on  the 

estimated  value.  However  v(i,j)  has  a  covariance  matrix  of 
2  2 

order  (N  xN  )  which  is  complicated  and  its  knowledge  is  not 


assumed  in  the  forthcoming  derivation.  Its  effect  on  the 
estimator  however  will  be ■ checked  by  simulation. 

The  general,  "weighted"  least  squared  cost  function  is 
defined  as: 


N-l  N-l 

t  52  w..  ( A(j> . . 

i  =  0  j  =  0 


‘♦ij.e5' 


(III-2) 


The  weights  can  be  used  to  account  for  different 

A 

quality  of  the  measurements  ,  and  must  have  two  main 

features : 


W.  .  =  W. . 
13  3i 


W..  =  0  for  i=j 
J 


Minimizing  the  cost  function  J  with  respect  to  0  will 


give : 


N-l  N-l 


o  =  4s-  =  L  £  W..  <A$..  -  A<f>.  .  )C-2 - Jiil)  ( III-3 ) 

d0  j_Q  1]  1]  13  00 


This  is  as  far  as  one  can  go  without  specifying  the  trace 
function  explicitly. 

To  solve  the  minimization  problem  two  approaches  are 


possible: 


-  For  "analytic"  trace  functions,  to  derive  an  explicit 


estimator. 
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-  For  "random”  trace  function  to  conduct  an  exhaustive 
or  "educated"  search  (which  is  not  efficient) ,  or  use 
some  form  of  gradient  search  algorithm  to  achieve 
minimization . 

In  this  dissertation  the  first  approach  will  be  followed 
due  to  the  fact  that  the  circular  and  linear  arrays  have 
"analytic"  trace  functions. 

1 .  Circular  Arrays 

The  trace  function  for  circular  array  is  (Eq.  II-8): 


(III-4) 


where  K  is  a  constant  for  a  given  trace  function. 
Inserting  the  above  equation  in  (III-2)  and  (III-3)  gives: 


J  = 


N-l  N-l 

£  £  wij  {Mij 


i=0  i=0 


K  sin[^-(i+j  )-0]  sin[j|(  i-  j  )  ] } 2 


(III-5 ) 


and 


N-l  N-l 


=  2K  2  1C  - 


i=0  j=0 


H  il 


K  sin  C^(i+j)-0]  sint^(i-j) 


x  sint^(i-j)]  cos  [j|(  i+ j  ) -0  ]  (III-6) 
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0  =  2K 


N-l  N-l 

E  E  j A4> ^ -  sin  ) 3cosC^(i+j  )  -  8] 

i  =  0  j  =  0 

-  KW^sin[jy(i+j  )-0]cos[g(i+j  )-0]sin2[jjy(i-j  )  ] 


N-l  N-l 

*  E  EWiji*ij<sil'Ti 

i=0  j=0 


•  2  7T  .  v  . 

smj]  )cos9 


N-l  N-l 


W.  -A4> .  .  (cos-^-i 
13  13  N 


2  IT 

cos^y-j  )  sin8 


K  v ^  /  •  4rr .  .  kn.  »  .  2n.  2tt 

-  £  >  >  wij<siniri  +  smT3-2sinTi  cos-p 


^  La  La 

i=0  j  =0 


-2cos  lyLi  sin-^-j  )cos20 


N-l  N-l 

SWij(coslfi  +  cosif^'2005!?1  cosTTj 

i=0  j=0 


+  2sin -^pi  sin-ypj  )sin20 


( III-7 ) 


In  general  the  last  two  terms  of  the  equation  are  non-zero; 
however,  for  some  important  and  practical  choices  of  they  are 
zero. 

Two  cases  for  which  the  last  two  terms  are  zero  are: 

-for  W. .  =1  for  all  i*j 
13 
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V 


-for  Wj'.  being  equal  in  groups  suuh.  that  each  pair  which 

2  IT 

has  the  same  value  for  cos  [-^-(i- j )  ]  will  have  the 

same  W.  .  . 

13 

The  importance  of  the  second  case  is  due  to  the  fact  that 
pairs  which  are  at  equal  distance  in  physical  space  will  have 
the  same  weights. 

For  the  cases  when  the  last  two  terms  of  equation  (III-7) 
are  zero  the  estimator  will  have  a  particularly  simple  form: 


A 


N-l  N-l 

£  £ 

tg~i 


^  O  *rr 

W.  .  A<J> .  .  (sin4r-i 
i]  N 


.  2tt  .  . 
sin^) 


N-l  N-l 

£  £ 

i=0  j  =  0 


„  ,  2tt. 

W.  .  .  .  (cOS-rfl 

13  13  N 


2rr .  . 
cos-jj-]  ) 


(III-8) 


In  the  general  case  however,  when  the  last  two  terms  are 
not  zero,  the  derivation  is  different.  One  can  write  the 
following  equation  from  (III  —  7) : 

0  =  cosQ  -  K2  sin9  -  cos29  +  sin20  CIII-9) 

where  K^,  K2 ,  K3,  are  expressions  which  can  be  easily 
calculated  for  a  given  estimated  trace  function. 


After  some  manipulations  one  obtains: 


which  is  an  equation  in  one  variable  and  can  be  written 

as: 

0  =  K±/  1+x2  -  K2x/l+x*  -  K3(1-x2)  +  2K4x  (III-ll) 

and 

0C=  tg-1  x  ( III-12 ) 

The  above  equation  can  be  solved  by  numerical  methods.  As 
mentioned  before  this  procedure  is  required  only  if  the  weights 
are  almost  arbitrary. 

An  additional  feature  of  the  estimator  in  CIII-8)  is  that 
it  does  not  require  the  knowledge  of  the  constant  K,  which 
contains  the  frequency  F,  diameter  D  and  speed  of  sound  C. 

K  =  2ttFD/C 

In  fact  it  is  important  to  realize  that  after  estimating  0 , 

X  can  be  estimated  and  as  a  result  the  speed  of  sound  C,  can 
be  calculated  when  the  frequency  F  and  the  diameter  D  are 
known.  This  feature  is  important  for  some  applications  such 
as  seismic  signal  processing  where  there  exists  the  problem 
of  dispersive  propagation  in  which  the  speed  of  sound  is 
trying  with  frequency.  The  estimation  of  K  will  be  done  as 


using  the  same  cost  function  (Eq.  III-5)  but  minimizing  with 

A 

respect  to  K  when  0  is  known: 


N-l  N-l 


:3J  _ 


{AA . .-K  sin[£(i+j  ) -0 ]sin[«-( i- j  )]} 


IT,  • 


2 


■N 


N 


i=0  j  =  0 


x  sin[^-(  i+j  )-0  ]sin[|j-(  i-j  )  ]  (III-13) 


and  the  estimator  is: 


N-l  N-l 


K 


y]  yT"]  sin[gCi+j)-9]  sintj(i-j)] 

i=0  j  =0 

Nri~ta - - - 

^  ^  ^  sin2[^(i+ j  )-0]sin2[^-(i-j  )  ] 

i=Q  j=0 


( 111-14 ) 


knowing  K: 


C  = 


2ttFD 

A 

K 


( III-15 ) 


2 .  Linear  Arrays 

Following  the  same  development  as  used  for  the  circular 
array,  the  trace  function  is  (Eq.  11-10): 


A*iij,e  =  Hp  e 
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=  KL(i-j)sin0  (III-16) 

N-l  N-l 

J  =  [A4>ij  -  KL(i-j)sin9]2  (III-17) 

i=0  j  =  0 


and 


0  =  IF  =  2Kl 


N-l  >1-1 

E  E 

i=0  j=Q 


W  .  .  C A4> .  .  -K.  (i- j  )sin8  ]  (i-  j  )cos6 
13  13  Li 


( III-18 ) 


and  N-l  N-l 

E  E'w1-” 

\  *  sln'1  - jpi - -  (III-19) 

n.  E  Ev^2 

i=0  j  =  0 


3 .  The  Degenerate  Trace  Function  Estimator 

The  degenerate  trace  function  for  circular  arrays  is 
(Eq.  11-13): 

A4>\  .  a  =  cos  (^-i-8)  CIII-20) 

13,0  C  N 

where  it  was  assumed 
* 

j  =  (i+  N/2 )  mod  N  (III-21) 
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Inserting  this  equation  to  the  estimator  (Eq.  III-8) 


and  assuming  all  equal  to  unity: 

N-l 

^  1  sin-^i 

«  .  tr-l  i=0 

eC.D.  ‘  tg  Fm - ^ -  (III-22) 

^  A4>£  cos-^-i 

i=0 


which  is  the  phase  of  the  first  harmonic  of  the  discrete 
fourier  transform  of  A$^. 


C .  CHAPTER  SUMMARY 


In  this  chapter  the  MMSE  estimator  of  the  bearing 
was  developed  for  circular  and  linear  arrays  based  on  the  trace 
function  principle.  Two  general  remarks  should  be  mentioned: 

1.  The  MMSE  estimator  is  equivalent  to  the  Maximum 
Likelihood  estimator  when  the  noise  is  gaussian  and  indepen¬ 
dent  among  the  data  points.  This  case  can  be  considered  the 
"best"  case  and  it  will  be  used  to  develop  the  Cramer-Rao 
Lower  Bound  for  the  estimator.  This  CRLB  will  serve  as 
a  benchmark  to  the  performance  of  the  estimator  which  was 
derived  by  simulations. 

2.  If  all  the  data  points  in  a  trace  function  are  esti¬ 
mated  from  the  same  time  frame/ s  the  two  halves  of  the  trace 


function  are  totally  dependent,  because 


A $ .  .  =-A$  .  .  result 
13  3i 


from  the  same 'heasurement" .  In  addition  the  terms  A$k^  are 
not  contributing  to  the  bearing  estimation  process. 
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As  a  result  of  these  the  number  of  contributing  terms  in 
the  trace  function  data  matrix  is: 


(N2  -  N)/2. 


The  summation  limits  in  the  estimators  must  be  changed 
accordingly: 

Equation  (III-8)  becomes: 


N-2  N-l 

E  E 

i=0  3=1+1 


.  .  ,  .  2n  •  .  2ir .  N 

W .  .  Ad> .  .  ( sxn-n-i-sm-rr-i ) 
lj  rij  N  N  J 


~  ,  .  2tr .  .  2ir .  v 

W  -  .  A«b  -  .  (sin-rT-i-sin-rr-n  ) 
ij  ri]  N  N  J 


(III-23) 


Equation  (III-19)  becomes 


N-2  N-l 

E  E 


2  -  gin"1  1=0  .  j  = 

0L  S  N-2  N-l 


WijA<J>ij(i-j) 


E  E 

i=0  j  =i+l 


W±j  (i-j)' 


( III-24 ) 


Equation  (III-22)  becomes 

N/2-1 


2  -  fa"1  1=0 

9C.D.  ‘  tg  N/2-1 


^  4*i  “"TT1 


E 

i=0 


*2  2 7T  . 

cos-jj-i 


(III-25) 
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However,  if  the  measurements  for  A<{) .  .  an  A<1>..  are  done 

13  31 

from  independent  data  blocks  the  equations  should  remain  in 
the  original  form. 

The  exact  structure  of  the  weights  was  not  treated 
here  and  it  will  not  be  pursued  in  this  work,  however,  it 
should  be  emphasized  that  this  is  a  subject  for  future  work. 
Especially  it  is  recommended  to  check  the  derivation  of  the 
weights  with  respect  to  their  connection  to  the  coherence  of 
the  two  channels  participating  in  each  pair. 

In  the  next  chapter  the  performance  of  the  estimators 
which  have  been  derived  in  this  chapter  will  be  checked  by 
simulation  to  demonstrate  the  effect  of  changes  in  various 
significant  parameters  on  the  estimator  performance. 


86 


IV.  ESTIMATOR  PERFORMANCE 


A.  INTRODUCTION 

In  the  previous  chapter  bearing  estimators  for  circular 
and  linear  arrays  were  developed  using  trace  function  prin¬ 
ciples  . 

In  this  chapter  the  performance  of  the  estimators  will 
be  checked  by  simulations .  The  two  main  measures  of  perfor¬ 
mance  are : 

-  Bias 

-  "Beamwidth"  which  here  is  defined  as  twice  the  standard 
deviation  of  the  bearing  estimate  (2 cj) 

In  addition  to  the  above  a  phenomenon  called  threshold 
effect  was  observed  which  is  expected  to  be  present  in  most 
nonlinear  phase  estimating  processors  [38,  p.  287].  This 
threshold  was  also  checked  for  various  parameters,  and  it  was 
defined  arbitrarily  as  the  input  phase  difference  variance 
level  for  which  the  output  "beamwidth"  is  less  than  five 
degrees . 

The  input  to  the  estimator  was  the  "estimated"  phase 
difference  .trace  function  (Eq.  III-l): 


1 


Two  noise  models  were  used  for  v(i,j): 

-  independent  gaussian  noise  from  pair  to  pair 

-  correlated  gaussian  noise  from  pair  to  pair 

The  model  of  the  correlated  noise  is  defined  in  Appendix  C 
along  with  the  simulation  method  and  computer  program.  The 
intent  of  including  the  correlated  noise  was  to  give  a 
feeling  about  its  effects  on  the  performance  and  not  so  much 
on  its  accurate  modeling. 

The  simulation  runs  were  performed  mainly  for  circular 
arrays,  however  some  performance  measures  were  checked  for 
linear  arrays  as  well  as  for  comparison  purposes. 

For  each  parameter  change  48  computer  runs  were  done  in 
order  to  get  reasonable  statistics.  The  mean  and  the  variance 
of  the  results  of  the  48  runs  were  calculated  and  they  served 
as  the  estimated  bearing  and  bearing  variance. 

B.  CRAMER-RAO  LOWER  BOUND  FOR  THE  IDEALIZED  ESTIMATOR 

As  a  benchmark  for  the  performance  of  the  estimators  the 
Cramer-Rao  Lower  Bound  (CRLB)  was  used.  It  was  calculated 
for  the  maximum  likelihood  estimator  which  is  equivalent  to 
the  MMSE  estimator  for  uncorrelated  gaussian  noise 
[23,  p.  198].  This  case  can  be  considered  as  "best  case"  and 
it  can  be  used  as  a  lower  bound  for  any  realizable  estimator. 
The  CRLB  for  a  maximum  likelihood  estimation  of  a  parameter 
"a"  is  (38;  p.  275,  41,  p.  335). 
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-  \ 


where  No  -  the  input  noise  level  (variance) 

s(t,a)-  the  "clean"  signal 

a  -  the  parameter  to  be  estimated 

Adaptation  of  the  above  to  the  present  case  will  be: 

2 

-  No  *  0,2  -  variance  of  the  estimated  phase  difference 

4*ij 

-  S(t,ot)  is  the  discrete  trace  function  fTr(i,j,0) 

-  The  integration  becomes  double  summation 
and  for  a  circular  array: 


°9e  -  - FFT? - KPT 


2 

a  * 

A*n 


2  ^2  X]  {^{K  sinc5<i+j)'e]sini:S'(i“j)3}} 

i=0  J=i+1 


JA<J> 


tPT  jjti- 

2K2  XI  XI  cos2[|-(i+j  )-0]sin2[^-(i+j )] 

i=  0  j-i+1 


(IV-3) 


for  a  linear  array: 


9L  ~ 


±1 


FP1  FPT 

2  1C  iC  {-^CKLsin8(i-j)]} 
i=0  j=i+l 

2 


A$ 


13 


£ 

i=0  j  =  i+l 


(IV-4) 


„  20  ->2 
2K^cos  0 
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The  above  bounds  were  plotted  along  with  the  simulation 
results  for  comparison  purposes  and,  as  it  will  be  shown  the 
simulation  results  achieve  this  bound  asymptotically. 

C.  BIAS 

1,  Circular  Array 

For  this  performance  test  the  simulation  was  run  for 
4,  8  and  16  element  5  meter  diameter  circular  arrays  at  100  Hz 
and  target  bearing  30°.  For  each  element  number  change,  the 

2 

input  phase  difference  variance  was  changed  between  4.35Rad 
(120°St  .Dev .  'j  down  to  6xl0-^  Rad^  (  .  45°St  .Dev  . )  in  25  steps. 

The  bias  was  calculated  for  each  run  by  the  equation: 

/N  A 

b  =  e  -  e 

A 

where  B  -  bias 

A 

0  -  estimated  bearing 

0  -  real  bearing  =  30° 

The  results  are  plotted  in  Figure  IV-1. 

The  "lock-on"  range  was  marked  on  the  plots  and  although 
the  bias  was  not  the  criterion  for  the  "lock-on"  range  defini¬ 
tion,  it  is  clearly  shown  that  the  bias  fluctuation  becomes 
very  small  within  the  "lock-on"  range. 

Observations : 

a.  Within  the  "lock-on"  range  the  bias  is  very  small 
and  is  bipolar,  a  fact  which  suggests  that  the  apparent  bias 
is  due  to  an  insufficient  number  of  runs  for  ehsemble  averaging 
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and  that  the  estimator  itself  is  really  unbiased . 

b.  Even  if  "a"  is  not  true  the  estimator  is  at  least 
asymptotically  unbiased. 

c.  The  correlated  measurement  noise  has  no  apparent 
effect  on  the  bias  compared  to  the  uncorrelated  noise. 

d.  As  the  number  of  elements  increases  the  fluctuations 
in  the  bias  are  smaller. 

2 .  Linear  Array 

For  this  performance  test  the  simulation  was  run  for 
a  linear  array  of  8  elements  with  1.91  meter  element-spacing 
at  100Hz  (taking  a  5  meter  diameter  circular  array  with  the 
same  number  of  elements  and  opening  it  up).  However,  this 
time  the  bearing  was  varied  as  a  parameter  because  unlike  the 
circular  array  the  linear  array  is  not  symmetric  in  the  x-y 
plane.  The  range  of  the  input  was  the  same  as  for  circular 
array.  Figure  IV- 2  shows  the  resultant  plots. 

In  addition  to  the  observations  of  the  circular  array,  it 
can  be  seen  that  the  bias  fluctuations  increase  toward  end- 
fire  but  even  at  those  bearings  the  bias  is  asymptotically  zero. 

3 .  Conclusions  With  Respect  To  Bias 

a.  It  can  be  assumed  that  the  MMSE  -  estimators 
developed  in  this  work  are  unbiased  or  at  least  asymptotically 
unbiased . 

* 

b.  The  lock-on  range  is  a  good  measure  for  the  region 
where  the  bias  fluctuations  settle  down. 

c.  Increasing  the  number  of  elements  reduces  the 
bias  fluctuations. 
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D.  BEAMWIDTH 


The  beamwidth,  which  was  defined  as  twice  the  standard 
deviation  of  the  bearing  estimates  (20q),  is  by  far  the  most 
important  performance  measure.  As  a  result  it  was  most 
extensively  checked  both  by  varying  the  various  parameters 
which  effect  the  performance  and  comparing  to  the  CRLB.  In 
this  section  most  of  the  simulation  runs  were  performed  for 
both  correlated  and  uncorrelated  "measurement"  noise  along 
with  the  CRLB. 

The  model  used  for  the  correlated  measurement  noise  is 
presented  in  Appendix  C.  This  case  of  correlated  noise  is 
the  actual  one  when  the  phase  difference  estimation  for  all 
pairs  is  done  simultaneously  (from  the  same  time  frame). 
However,  if  the  estimation  will  be  done  from  different  time 
frames  for  various  groups  of  pairs  the  measurement  noise  can 
be  accounted  as  uncorrelated  and  then  the  uncorrelated  mea¬ 
surement  noise  case  will  apply.  The  effects  on  the  perfor¬ 
mance  of  both  cases  is  shown  throughout  this  section. 

1 .  Circular  Arrays 

a.  Beamwidth  Versus  Phase  Difference  Estimation  Noise 
Variance 

Figures  IV-3  through  IV-6  show  the  results  of  the 
simulation  runs  for  this  performance  measure  for  4,  8,  12  and 
16  element  circular  arrays  of  equal  diameter  of  5  meters  at 
100Hz.  For  the  4  and  8  element  arrays  both  correlated  and 
uncorrelated  measurement  noise  is  given.  For  all  the  plots 
the  CRLB  is  shown  along  with  the  simulation  results  as  a 
benchmark  to  performance. 
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Observations : 

-  In  all  plots  the  results  and  the  CRLB  exhibit  a 
clear  threshold  effect. 

-  The  correlated  noise  assumption  (which  is  the 
realistic  one)  does  not  have  a  major  effect  on  the 
performance . 

-  Under  both  correlated  and  the  uncorrelated  measure¬ 
ment  noise  assumption  the  estimator  achieves 
asymptotically  the  CRLB  at  low  level  input  noise 
variance,  thus  the  estimator  is  asymptotically 
efficient  [38 ,  p.  276]  for  decreasing  input  noise 
level. 

b.  Beamwidth  Versus  Numbers  Of  Elements  At  Constant 
Diameter 

Figure  IV-7  shows  the  dependence  of  the  beamwidth 
on  the  number  of  the  elements  for  a  circular  array  of  5  meter 
diameter  at  100Hz.  The  plots  were  done  for  three  representa¬ 
tive  input  variance  levels.  Only  the  uncorrelated  measure¬ 
ments  along  with  the  CRLB  are  shown  which  are  adequate  to 
represent  the  basic  behavior. 

Observations : 

-  A  saturation  effect  is  clear  for  all  three  variance 


levels  shown;  at  each  noise  level  at  a  given  diameter 
and  frequency  there  is  a  limit  on  the  number  of 
elements  required,  increasing  the  number  of  elements 
above  this  "saturation  number"  will  improve  perfor¬ 


mance  only  marginally.  For  example  for  the  lowest 

-3  2 

curves  (a."?  =2.25x10  Rad  )  this  number  is 

AO  .  . 

13 


99 


approximately  16  elements,  for  higher  noise  levels 
it  will  be  higher. 

-  Again  the  simulated  estimator  performance  approaches 
the  CRLB  asymptotically  for  large  number  of  elements. 

c.  Beamwidth  Versus  Diameter  For  A  Constant  Number 
Of  Elements 

Figure  IV-8  shows  the  dependence  of  the  beamwidth 
on  the  diameter  of  a  circular  array  of  8  elements  at  100Hz. 

The  plots  were  done  for  two  representative  input  variance 
levels.  The  correlated  and  uncorrelated  measurements  along 
with  CRLB  are  shown. 

Observations : 

-  A  saturation  effect  is  clear  in  this  case  too; 
increasing  the  diameter  of  the  array  above  certain 
limits  will  improve  performance  (smaller  beamwidth) 
only  marginally.  This  limit  appears  to  be  approxi¬ 
mately  at  15  meter  diameter  which  is  one  wavelength 
at  100Hz. 

-  The  behavior  of  the  simulation  resultant  curves  is 
similar  to  the  CRLB  and  they  approach  the  CRLB 
asymptotically . 

2 .  Linear  Arrays 

The  analysis  of  performance  for  the  linear  array  is 
much  less  extensive  than  for  the  circular  array  and  its  aim 
is  two  fold: 

-  It  is  the  most  standard  geometrical  configuration 
treated  in  literature  and  a  reader  of  this  work 
can  readily  compare  it  with  other  results . 
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-  It  was  compared  to  the  circular  array  in  some 
aspects  of  its  performance. 

In  this  section  only  beamwidth  versus  phase  difference 
estimation  noise  variance  plots  are  presented.  Figures  IV-9 
and  IV-1Q  show  the  resultant  curves  for  a  linear  array  of 
5  meter  length  (same  aperture  as  the  circular  array)  with 
4  and  8  elements  at  100Hz.  The  plots  are  for  both  correlated 
and  uncorrelated  measurements  along  with  the  CRLB. 

Figure  IV-11  and  IV-12  show  the  results  for  a  linear 
array  of  8  elements  and  1.91  meter  interelement  spacing. 

This  configuration  is  equivalent  to  an  "opened  up"  version 
of  a  5  meter  diameter  circular  array  preserving  the  inter¬ 
element  distance  and  number  of  elements.  Figures  IV-11  and 
IV-12  differ  in  the  look  direction  which  causes  a  different 
performance,  as  expected,  because  the  asymmetry  of  the  linear 
array  in  the  x-y  plane . 

Observations : 

-  The  linear  array  estimates  present  a  similar  be¬ 
havior  as  the  circular  array  with  respect  to  their 
performance  relative  to  the  CRLB  -  they  are 
asymptotically  efficient. 

-  A  linear  array  with  the  same  aperture  and  number 
of  elements  as  a  circular  array  has  a  slightly 
worse  performance,  the  reason  being  the  trend  to 
concentration  of  elements  at  the  ends  for  the 
circular  array,  a  characteristic  which  improves 
bearing  estimation. 


Beamw^dth 


Figure  IV- 


Beamwidth 


Beamwidth 


Input  variance  Figure  IV-11 


A  linear  array  with  the  same  number  of  elements 

and  linear-element  distance  has  however  a  much 

better  performance  than  the  equivalent  circular 

array  because  of  the  much  bigger  aperture. 

For  a  linear  array  the  beamwidth  increases  toward 

endf ire ,  a  fact  which  is  expected  from  the  equation 

2 

of  CRLB  (Eq.  IV-4)  where  the  dependence  on  1/cos  9 
is  observed.  The  simulation  results  show  the  same 
behavior . 
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E.  THRESHOLD  EFFECT  AND  ITS  IMPORTANCE 

As  mentioned  in  the  introduction  to  this  chapter  a 
threshold  effect  which  can  be  expected  in  this  kind  of  non¬ 
linear  estimator  was  observed.  In  this  work  the  threshold 
was  defined  as  "lock-on  point"  and  it  is  the  value  of  the 
input  variance  for  which  the  output  beamwidth  becomes  less 
than  five  degrees.  The  number  of  five  degrees  was  chosen 
arbitrarily  and  it  can  be  any  number  in  the  vicinity  of  the 
knee  of  the  performance  curve  (see  for  example  Figures  IV-3 
through  IV-6  for  circular  arrays  and  Figure  IV-9  through 
IV-12  for  linear  arrays). 

The  range  of  input  variances  for  which  the  beamwidth  is 
less  than  the  threshold  value,  in  this  case  5  degrees,  is 
defined  as  "lock-on  range". 

The  importance  of  this  issue  is  due  to  the  following: 

-  The  lock-on  point  can  be  translated  to  a  threshold  on 
the  coherence  at  a  specific  frequency  and  only  those 
phase  difference  measurements  for  which  this  threshold 
is  exceeded  are  processed  by  the  trace  function  estima¬ 
tor.  (This  connection  will  be  further  explained  in 
Chapter  V.)  This  thresholding  may  ensure  that  only 


those  frequencies  which  will  guarantee  a  certain  beam- 
width  will  be  processed.  It  should  be  mentioned  that 
the  threshold  is  frequency  dependent  exactly  in  the 
same  manner  as  the  diameter  dependence  which  will  be 
shown  later  in  this  section. 

-  As  it  can  be  seen  from  Figures  IV-1  and  IV-2  operating 
within  the  "lock-on  range"  ensures  also  that  within 
the  beam  width  there  will  be  less  variability  of  the 
resultant  bearing  estimate. 

The  dependence  of  the  lock-on  point  for  a  circular  array 
on  number  of  elements  and  diameter  is  shown  in  Figures  IV-13 
and  IV-14  respectively.  The  plots  are  for  uncorrelated 
measurements  and  for  the  CRLB. 

The  frequency  dependence  of  the  lock-on  point  is  exactly 

the  same  as  for  diameter  and  is  not  shown  separately  (both 
2 

have  an  "x  "  effect). 

The  values  of  SNR  shown  on  the  plots  are  obtained  from 
variance  values  related  back  to  the  system  input  and  their 
aim  is  to  give  a  feeling  of  the  change  in  the  noise  tolerance 
of  the  system  as  the  diameter  (or  frequency)  and  number  of 
elements  are  changed. 

Although  the  simulation  results  do  not  achieve  CRLB 
within  the  range  of  values  presented,  nevertheless  the  trend 
shown  suggest  that  it  is  eventually  asymptotically  achieved 
for  higher  number  of  element  and  larger  diameter. 
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Figure  IV-14 


F.  PERFORMANCE  SUMMARY 


1.  Within  the  "lo<'k-on  range"  the  estimator  is  practically 
unbiased. 

2.  With  respect  to  the  beamwidth  performance  the  estima¬ 
tor  is  asymptotically  efficient  i.e.  it  asymptotically  achieves 
the  CRLB  when  the  input  variance  decreases,  or  the  number  of 
elements,  diameter,  or  frequency  increases. 

3.  A  threshold  effect,  defined  here  as  "lock-on"  was 
observed  and  it  enables  one  to  tailor  the  desired  system  per¬ 
formance  . 

4.  There  is  a  limit  on  the  number  of  elements  at  a  fixed 
diameter  (or  on  the  diameter  at  a  fixed  number  of  elements) 
beyond  which  the  beamwidth  will  only  improve  marginally  with 
the  increase  of  number  of  elements  or  array  diameter  respec¬ 
tively. 

5.  The  trace  function  estimator  applies  to  circular 
arrays  as  well  as  to  linear  arrays ,  however  for  the  same 
aperture  and  number  of  elements  the  circular  array  performs 
slightly  better. 

In  the  next  chapter  on  system  considerations  a  system 
which  uses  the  estimators  developed  here  will  be  suggested. 

The  problems  involved  in  implementing  such  a  system  will  be 
analyzed. 
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V.  SYSTEM  CONSIDERATIONS 


A.  INTRODUCTION 

The  motivation  behind  this  work  has  been  to  find  a  solu¬ 
tion  to  the  practical  problem  of  improved  bearing  estimation 
using  a  small  array.  The  trace  function  concept  offers  a 
way  to  solve  the  problem  but  is  not  the  solution  itself.  In 
order  to  solve  the  problem  giving  a  meaningful  and  practical 
result  the  "pieces"  must  be  connected  together  into  a  system 
concept  which  considers  the  signal  from  the  sea  to  the  display. 
In  this  chapter  such  a  possible  system  will  be  suggested  and 
outlined.  It  is  not  the  intent  to  design  a  system  but  rather 
to  highlight  those  points  which  are  to  be  considered  in  a 
system  using  the  trace  function  concept. 

The  main  task  would  be  to  process  the  received  data  at 
the  sensors  (hydrophones)  to  bring  it  to  a  trace  function 
format . 

Once  the  data  is  in  phase  difference  and/or  time  delay 
trace  function  format  the  estimation  process  developed  in 
previous  chapters  can  be  applied  and  the  bearing  estimated 
for  each  frequency  resolution  cell  of  interest  or  for  the 
wide  band  signal. 

In  the  process  of  performing  this  task  there  is  a  major 
difficulty  to  be  overcome;  any  processing  scheme  used  must 
preserve  the  relative  phase  between  the  signals  received  the 
same  as  they  were  at  the  receiving  sensors  in  order  to  result 


in  a  meaningful  trace  function. 
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Two  process  steps  may  be  identified  which  pose  this 
difficulty: 

-  The  decorrelation  (prewhitening)  both  in  space  and  time. 

-  The  phase  unwrapping  which  is  needed  because  of  the 
inherent  2ir  ambiguity  present  in  any  phase  measuring 
system. 

Those  and  other  points  will  be  reviewed  and  analyzed  briefly 
in  this  chapter.  Some  suggestions  will  be  made  to  overcome 
the  mentioned  difficulties. 

In  Figure  V-l  a  conceptual  system  functional  block  dia¬ 
gram  is  presented.  It  applies  the  trace  function  principles 
as  discussed  in  previous  chapters.  In  the  forthcoming  sections 
the  functions  and  problems  associated  with  various  blocks  will 
be  treated . 

B.  SENSOR  ARRAY 

The  sensor  which  will  usually  be  a  hydrophone  is  the  first 
part  of  the  system  interacting  with  the  signal.  The  major 
characteristics  which  have  to  be  considered  at  the  sensor 
array  are: 

-  Accurate  and  known  geometrical  configuration. 

-  Equal  phase  response  of  the  sensors  across  the  frequency 
band  of  interest. 

If  the  phase  response  characteristic  cannot  be  satisfied 
because  of  practical  problems,  at  least  an  accurate  response 
chart  at  all  frequencies  of  interest  must  be  obtained  and 
stored  in  a  computer  memory.  This  data  can  be  used  in  turn 
to  compensate  for  the  different  response  characteristics  of 
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the  various  sensors  when  the  phase  difference  calculations 
are  performed  for  each  pair  of  sensors.  This  way  the  rela¬ 
tive  phase  differences  at  the  input  to  the  system  can  be 
preserved.  Any  other  change  in  response  of  the  sensors  will 
be  treated  by  the  system  as  additional  noise  on  the  phase 
difference  measurements. 

C.  PREPROCESSING 

In  this  part  of  the  system  preamplification,  bandlimiting 
and  A  to  D  conversion  takes  place.  The  only  requirements  are: 

-  zero  phase  or  at  least  linear  phase  response  at  all 
frequencies  of  interest. 

-  Equal  phase  response  for  all  channels. 

After  the  A  to  D  conversion  digital  methods  will  be  used 
and  the  preservation  of  equal  phase  response  becomes  much 
simpler. 

D.  SPATIAL  DECORRELATION  (WHITENING) 

The  processing  of  the  signals  received  by  an  array  of 
sensors  in  a  noisy  background  is  based  upon  the  relative 
coherency  of  the  directional  signal  compared  to  the  relative 
incoherency  of  the  background  (ambient)  noise.  When  this 
difference  is  distinctive  (i.e.  the  noise  is  almost  uncor¬ 
related  from  sensor  to  sensor)  the  processing  is  relatively 
easy  and  not  very  much  dependent  on  the  noise  model.  This 
characteristic  is  achieved  when  the  separation  between  sensors 
is  in  the  order  of  a  half-wavelength  for  the  frequencies  of 
interest.  The  task  is  very  much  complicated  where  the  separation 
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is  small,  because  the  high  correlation  of  the  background  noise 
narrows  the  difference  between  the  signal  and  noise  in  that 
domain. 

In  the  problem  defined  in  this  dissertation  the  case  is 
of  small  arrays  where  the  whole  array  aperture  is  much  smaller 
than  a  wavelength  at  the  lowest  frequencies  of  interest  and  as 
such  the  noise  will  be  highly  correlated. 

The  conventional  way  to  decorrelate  the  noise  among  the 
sensors  is  by  multiplication  with  the  inverse  noise  covariance 
matrix  or  the  inverse  signal  and  noise  covariance  matrix. 

This  operation  is  prescribed  by  every  "optimal  array  processor" 
(2,  5,  39)  and  it  can  be  performed  based  on  a  priori  know¬ 

ledge,  "estimate  and  plug",  or  adaptive  calculation  of  the  in¬ 
verse  matrix. 

This  operation  is  inherentely  not  phase  preserving  and 
this  fact  will  be  shown  in  the  forthcoming  simple  example: 

Assume  two  sensors  with  the  signal  part  of  the  input 
being : 


iwt 

X1 

_X2_ 

ei(wt+A(}>) 

clearly  the  phase  difference  is  A$. 

A  typical  covariance  matrix  has  the  form: 


Q  = 


1 

P 


P 

1 


( V  —  2  ) 
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then 

Q 

where 


-1 


1-P‘ 


1  -p 

and  Q_is  =  i 

1  -a 

-p  1 

D 

_-a  1_ 

(V-3) 


a  =  Iz^ZfL 
P 


and  i 


^a-p2)(i-A-p2) 

The  prewhitening  process  is  defined  as  Q  ^x  as  it  can  be  seen 
from: 


and 


q-1’;  -  4 

D 


;}=  Q"*5  E  {x  x11}  Q~h 

=  Q"35  Q  Q 

-alfe-  1. 

1  fe1"1 

-a  1 


i(uit+A4> ) 


h*eiwt+eiCwt+A*) 


iu>t  .  A. 

xx  =  (1  -  aelA*) 


„  l(i)t  .  .  , 

x2  =  £3-  (elA*-a) 


(  V-4 ) 

(  V-5 ) 


After  some  manipulations  the  phase  difference  between 
x2  and  x1  is: 

/ 

~  =  tg'1  sinA4>(l-a2)  _  -1  sin^  ( P  2-l+  A-o2] 

*2*1  cosA<Ji(l+a  )-2a  x  (cosA<J>-p)(l~  >4-p^) 

(V-6 ) 

which  clearly  cannot  be  equal  to  A$  for  any  value  of  p  except 
p-*-0  which  the  trivial  uncorrelated  case. 

In  order  to  preserve  phase  from  the  total  system  point  of 
view  one  must  correct  for  the  above  phase  distortion. 
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A  new  method  called  "phase  preserving  spatial  decorrela 
tion  (or  prewhitenning)"  by  the  author  is  proposed  and  out¬ 
lined  here.  The  notion  is  that  the  phase  distortion  caused 
by  the  decorralation  process  can  be  calculated  and  compen¬ 
sated  for,  after  the  cross  spectral  estimation  process,  pro¬ 
vided  that  one  knows  exactly  the  inverse  covariance  matrix. 

To  demonstrate  the  above,  the  previous  example  is  con¬ 
tinued  : 

the  original  phase  difference  was  A<f> 

xlx2 

the  resultant  phase  difference  after  decorrelation  was 


J-  -  =  tg"1  g4q.4i<Iz3 _ 

xlx2  cosA<j>(l+aZ)~2a 


The  difference  between  the  two  is  the  compensation  term 


A<J>-A<f>  =  tg 


-1  2asinA<j)(l-acosA({)) 


1-a  cos  2A4>-2acosA<f> 


( V-7 ) 


The  above  compensation  term  is  impractical  to  calculate  be¬ 
cause  it  is  dependent  on  the  real  phase  difference  which  is 
unknown.  As  a  result  no  simple  compensation  can  be  used. 
Equation  (V-6)  can  be  solved  numerically  however  for  A<j> 

xlx2 

once  A4>~  -  is  measured  provided  p  is  known. 
xlx2 

As  shown  by  this  simple  example  this  approach  is  very 
complicated  even  for  simple  arrays  (2  elements)  because  it 
has  to  be  calculated  for  each  frequency  bin  and  each  sensor 
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pair.  For  practical  number  of  elements  in  an  array  this  kind 
of  operation  becomes  prohibitive.  To  solve  the  problem  a 
different  approach  has  to  be  found. 

As  the  ambient  noise  is  a  wideband  process  its  cross 
covariance  function  has  a  limited  time  width 

Pij(T)  <  e  for  t  >  tq  (V-8) 

i.e.  it  can  be  made  arbitrarily  small  if  long  enough  delay 
between  the  elements  is  applied. 

Figure  V-2  shows  such  a  covariance  function  for  a  band¬ 
width  of  50  to  250  Hz  a  power  spectrum  of  -6dB/octave  for 
2.5  and  5  meter  distance  between  sensors. 

The  graph  clearly  shows  the  decay  of  the  covariance  "p" 
as  function  of  delay  time  "t".  This  characteristic  can  be 
used  for  the  decorrelation  process.  One  can  introduce 
successive  and  known  delays  between  the  elements  of  an  array 
such  that  the  smallest  delay  is  enough  to  decorrelate  the 
noise  down  to  a  desired  level  for  the  closest  pair  in  the 
array . 

After  this  process  the  signal  is: 


where  A  is  the  incremental  time  delay  introduced  between  the 
elements . 

It  is  obvious  that  to  compensate  for  this  process  the 
only  operation  to  be  done  is  to  re-add  at  each  frequency 
a  phase  compensation  term  to  the  calculated  phase  dif¬ 
ference.  This  term  will  be  for  the  pair  "i,j": 


comp. 


Mi-j  )«£ 


(V-10) 


This  method  will  not  affect  narrow  band  signals  but  may  also 
decorrelate  partially  the  wide  band  signal. 

This  effect  can  be  minimized  if  A  is  chosen  such  that  is 
large  enough  to  decorrelate  the  noise  but  (N-l)A  is  smaller 
than  the  decorrelation  time  of  the  wide  band  signal. 

In  order  to  perform  this  operation  one  must  have  the 
noise  and  signal  correlation  function  which  can  be  obtained 
either  from  apriori  knowledge  or  from  estimations  made  in 
situ  • 

The  above  method  should  be  viewed  as  a  suggestion  and  not 
a  firm  solution  mainly  because  it  is  based  on  theory  only 
without  simulation  or  measurement  back-up. 

E.  REVIEW  OF  SPECTRAL  ESTIMATION 

This  section  will  summarize  the  spectral  estimation  tools 
needed  for  subsequent  estimation  of  time  delay  and  phase 
difference. 
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As  it  will  be  shown  later  if  one  has  infinitely  long  data 
record  he  can  get  both  stable  estimates  and  infinitely  good 
resolution  of  the  power  spectra,  and  phase  spectra  of  a 
stationary  random  process.  The  main  factor  which  limits 
achieving  the  desired  performance  is  the  limited  length  of 
data  available.  The  reasons  are  many,  but  one  of  the  most 
important  ones  is  the  time  limit  on  the  stationarity  of  the 
environment.  In  the  case  of  data  in  the  ocean,  a  few  minutes 
can  be  considered  as  a  good  measure  for  stationary  environment. 

The  analysis  which  will  follow  reproduces  the  main  results 
in  Ref.  [25,  26,  27]  which  deal  with  the  limited  time  record 
auto  and  cross  spectral  estimations. 

The  two  basic  parameters  which  control  the  performance 
of  spectral  estimation  are  the  available  record  length,  T, 
over  which  the  sample  of  the  random  process  is  assumed 
stationary  and  the  desired  frequency  resolution,  B,  of  the 
spectral  analysis.  Large  values  of  BT  give  good  performance 
but  small  values  are  often  encountered  because  of  short  T 
available  or  requirement  for  fine  resolution  B. 

In  order  to  take  advantage  as  much  as  possible  of  the 
available  data  a  method  was  developed  [40]  to  obtain  good 
reduction  of  variance  by  overlapped  processing  of  the  data. 

This  method  is  less  efficient  than  the  same  number  of  segments 
non-overlapped  but  it  is  much  better  than  using  the  available 
data  in  non-overlapped  form. 

Several  questions  arise  when  one  attacks  the  problem: 
what  percentage  of  overlap  to  use?,  which  windows?,  how  big 
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the  FFTs  have  to  be?  Those  questions  and  others  were  answered 

in  Ref.  [25,  26]  and  the  results  follow. 

The  stationary  random  process  x  and  y  are  observed  for  a 

time  T  seconds,  thus  x(t),  y(t)  are  available  for  0<t<T.  One 

wishes  to  estimate  the  cross  spectrum  G  (f)  with  a  resolution 

xy 

of  B  Hz  over  a  bandwidth  of  ft  Hz. 

The  method  of  segmenting  the  available  data  is  shown  in 
Figure  V-3. 


Figure  V-3 

Segmenting  Data  With  Overlapp 

The  connection  between  the  various  terms  in  Figure  V-3  is: 

(P-1)S+L  _<  T  ( V-ll ) 

where  P  -  number  of  segments 

S  -  shift  of  each  adjacent  window 
L  -  data  window  length 
T  -  total  available  data  length  (time) 
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c.  Estimate  of  cross  spectra 


Gxy(fn> 


PTT  Z  W  VV 


vv°lpw 


(  V-16 ) 


( Im[G  (f  )]i 

Phase  Ph  ( f  )=tg_i< - =  **  (f  ) 

xy  n  / Re[G  (f  )])  xy  n 

xy  n 


( V-17 ) 


Amplitude  A  (f  )  =  /Re[G  (f  )])*+(Im[G  (f  )])  (V-18) 

*  xy  n  xy  n  xy  n 


d.  Estimate  of  magnitude  squared  coherence 


IVV1 


| G  (f  ) 
1  xy  n 


G  (f  )G  (f  ) 
xx  n  yy  n 


( V-19 ) 


NOTE:  The  magnitude  squared  coherence  function  whose 

role  will  be  discussed  later,  is  given  here  because 

of  its  connection  to  the  other  estimation  processes, 

2 .  Statistics  Of  The  Estimates 

The  spectral  estimates  presented  are  random  variables. 

Their  means  and  variances  will  be  stated  under  the  assumptions 

that  x  and  y  are  Gaussian  random  processes  and  the  frequency 

2 

resolution  of  the  spectral  window  |W|  is  narrower  than  the 
finest  detail  in  spectrum  G. 


where 


!n>  =  Z  M<3>  e'i2’,jn''L 


(V-20) 


127 


a .  Means 


E  <«xx'V>  3  Gxx(fn>  /  lW<v 


)  I  2dv 


(V-21) 


E  {G  (f  )>  =  G  (f 
xy  n  xy 


n)  /  jW(v)j2dv 


(V-22) 


One  can  assume  without  loss  of  generality  that  the  windows 
are  normalized  i.e.  J  |Wj2  =  1  then: 


E  {G  (f  )}  8  G  <f  ) 
xx  n  xx  n 


(V-23) 


E  {G  (f  )}  ~  G  (f  ) 
xy  n  xy  n 


( V-24 ) 


Thus  the  estimators  are  unbiased 


E  {A*xy(fn)}  =  AV(fn) 


(V-25) 


E  {Axy(fn)>  =  Axy  =  |Gxyj 


CV-26) 


E{|YV„|2}  =  1 Y  |  2  +  kl-|yi2)2(l+^L!)  ( V-27 ) 


xy 


(after  [28]) 


The  estimate  is  biased  by  the  second  term  of  the 
equation.  A  note  about  the  above  bias  is  given  below. 

It  is  calculated  for  nonoverlapping  segments.  However, 
when  overlap  occurs  this  bias  is  substantially  reduced  to  a 
value  of  about  50%  of  the  above  figure  at  50%  overlap. 
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b.  Equivalent  number  of  degrees  of  freedom 

The  convenient  way  to  describe  the  stability  of 
any  positive  estimate  is  its  equivalent  number  of  degrees  of 
freedom  (EDF)  [18].  The  definition  of  this  stability  constant 
is : 


edf  =  k  =  l[a-v.er*  se]2  =  ( V— 2  8 ) 

variance  Var{G(f)} 

Basically  using  tables  as  in  Ref.  [18]  one  can 
estimate  his  confidence  in  the  measurements. 

For  example  for  K=500  90%  of  the  measurements  will 
be  in  the  range  of  0.92  to  1.081  of  the  average.  Maximizing 
those  K's  will  be  one  of  the  important  criteria  for  the  over¬ 
lap  and  window  choice. 

c.  Variances 


(P-1) 

Var{Gxx<fn»  3  Gxx2(V  F  L  <l-4L>IVkS)|2  <V-29) 

k=-(P-l) 

« 

where 

*w<T>  =  J  w(t)w  (t-T)dt  (V-30) 

(P-1) 

toIVf,»  5  Gxy2(fn>F  Z  (l-4L>IVkS>|2  (V-31) 

k=-(P-l) 

At  this  point  we  can  define  EDF  (k)  for  both  auto  and  cross 
spectra  [25,  26]. 
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<t>  (kS ) 


for  auto  spectra  (V-32) 


E 


k=-(P-l) 


|  Y  |  K, 

'xy'  I 


for  cross  spectra 


(V-33) 


Var  {A  (f  )}  =  G  (f  )G  (f  )[l+|y  (f  )j2]/KA  (V-34) 

xy  n  xx  n  y^’  n  1  'xy  n  1  A 


1-|YXv(fn)l 
Var(A<D  (f  )}=  - %  - 

xy  n  |y  (f  ) | 2k. 

1 'xy  n  1  / 


( V-35 ) 


Those  are  based  on  the  assumption  that 

2 

K  SK.  |y„„|  »1*  It  should  be  noted  that  neither  of  the  two 
c  a  xy 

cross  variances  depend  on  the  actual  phase  and  the  phase 
variance  is  not  dependent  on  actual  data  values  at  all,  but 
only  on  the  coherence  between  the  two  sequences  and  the 
processing  scheme. 


Vart^xy(fn)|2>  5  FlTxy(fn)|2a'lYxy<fn)|2)2  for  V*0 


( V-36 ) 


3 .  Design  Tradeoffs 


As  one  can  see  from  all  the  above  statistics,  k  and  P 
are  playing  major  roles  in  the  variance  reduction. 

In  order  to  improve  performance  one  has  to  maximize 
k  and  P  (which  are  connected)  for  a  given  amount  of  available 
data. 

In  Ref.  [25]  this  subject  was  thoroughly  investigated 
for  various  window  functions. 

Following  is  a  summary  of  important  results: 

a.  The  EDF(k)  has  a  maximum  as  related  to  P  for  a 
given  amount  of  data  T. 

b.  This  maximum  was  proven  to  be  almost  independent 
of  the  window  employed  for  large  BT  products 

max  EDF  3  3(BT-1)  (V-37) 

c.  The  required  fractional  overlap  which  is  defined 
as 

FO  =  ^  =  -~|5^  ( V-38 ) 


is  significantly  increased  for  max  EDF  compared  to  0.99  max 
IDF  so  it  is  not  recommended  to  design  for  max  EDF,  but  1%  or 
2%  less. 

As  an  example  using  cosine  window  (Hanning)  and 
BT=64  for  max  EDF  P=146  segments  are  needed  but  for  0.99 
max  EDF  only  112  segments  are  needed.  Clearly  the  additional 
30%  of  processing  is  not  worth  the  one  percent  increase  in  EDF. 
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The  optimal  FO  comes  out  to  be  in  general  not  a 


convenient  number  to  use,  e.g..  0.61  for  cosine  window.  One 
will  probably  choose  the  closest  convenient  value  like 
F0=0.5  which  will  give  for  cosine  window  EDF=92%  of  max  EDF 
which  is  a  reasonable  figure. 

d.  As  the  max  EDF  is  practically  independent  of  the 
window  employed  no  tradeoffs  had  to  be  made  to  realize  better 
side  lobe  response  of  the  more  sophisticated  windows.  This 
issue  is  effecting  the  size  of  the  particular  FFT’s  to  be 
performed. 


N 


s 


(V-39) 


where 

N  =  number  of  samples  in  one  FFT 
s 

Cw  =  half  power  bandwidth  constant  of  a  window 
At  =  sampling  interval 

Clearly  for  a  "better"  window  Cw  is  bigger  in 

general  and  N  will  be  bigger.  Thus  the  side  lobe  performance 
s 

for  a  given  B  can  be  improved  by  paying  with  bigger  FFT's  and 

not  affecting  stability  of  the  estimates. 

4 .  Effect  Of  Pure  Tones  (Lines)  In  The  Spectrum 

The  results  shown  so  far  were  true  under  the  assumption 

2 

that  the  B  of  the  spectral  window  |W(f  )|  is  narrower  than 
the  finest  detail  of  the  spectrum.  For  tonal  signals  this 
will  not  be  the  case  and  a  change  in  the  statistics  will  occur 
but  the  processing  method  will  not  change.  In 


Ref.  [27]  those  changes  were  made  and  the  following  statistics 
derived : 

a.  Mean 


E{G(f  )} 
n 


G(f  )  .. 
n  wide 


G<fn>line 


(V-40) 


where 


G(f  )  ^  =  The  previous  wide  band  power  spectrum 


G(f  ). •  =  A2 I W ( f  -f  ) J  2 

n  line  1  no’ 


(V-41) 


A  -  line  amplitude 
fQ-  line  frequency 


E{G  (f  )}  =  G  (f  )  .  .  +A  A  |W(f  -f  ) 

xy  n  xy  n  wide  x  y1  no 


( V-42 ) 


b.  Variance 


Var  {G(f  ) }  =  Var  (G(f  ) }  .  .  +2G(f  )  . .  G(f  ) 

n  n  wide  n  wide  n  narrow 


P-1 


F  2  <1-LFi)1w<kS)ei 


2?rk(f  -f  )S 
n  o 


(V-43) 


k=l-P 


from  which  a  bigger  absolute  value  of  the  variance  can  be 

expected  when  a  tonal  is  present.  However,  if  the  EDF  (K) 

is  calculated  for  G<f  )  -j  <<G(f  )  line  which  is  the  case  in 

n  wide  n 

general  when  lines  are  present,  one  can  see  that  this  relative 
measure  is  increased  which  means  a  more  stable  estimation. 
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*1 


For  G(f  )  >>  G(f  )  .  .  : 

n  narrow  n  wide 


EDF 


^tonal 


PG(f  ) 

_ n  narrow _ 

SI,„U  T'  (l-1)i)*„<KS)ei2,,k(fn-fo,S 
k=l-P 

(V-44) 


this  ^  is  bigger  than  the  one  for  wide  band 

alone  by  a  factor  of  approximately 


%G(f  )  /G(f  )  .  .  >> 

n  narrow  n  wide 


( V-45 ) 


All  the  other  equations  for  cross  spectral  estimation  statistics 

will  still  hold  with  the  new  Ktona^* 

It  is  worthwhile  to  note  the  | y  \  will  be  much 

xy 

bigger  in  this  case  (very  close  to  unity)  and  its  bias  and  variance 
will  be  very  small  because  of  the  increased  "S/N". 

F.  PHASE  DIFFERENCE  ESTIMATION  AND  THE  UNWRAPPING  PROBLEM 

The  basic  phase  difference  estimation  is  done  entirely  by 
the  cross  spectral  estimation  process  as  reviewed  in  the  pre¬ 
vious  section.  However  there  are  some  additional  steps  to  be 
performed  in  order  to  construct  the  phase  difference  trace 
function  or  estimate  the  wide  band  time  delay  [32,  36]: 

-  The  frequency  bins  for  which  further  processing  will  be 
done  must  -be  detected. 

-  The  phase  difference  estimated  value  must  be  recompen¬ 
sated  for  the  decorrelation  process  as  stated  in  Eq.  (V-10). 


i 
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-  The  resultant  phase  difference  has  to  be  unwrapped 
because  it  is  measured  only  in  the  range  tv  while  it 
can  have  any  ±n2iT  ambiguity. 

1 .  Detection  Problem 

The  detection  of  the  frequency  bins  which  are  candi¬ 
dates  for  further  processing  is  important  for  the  following 
reasons : 


-  To  reduce  the  false  alarms  in  the  system. 

-  To  reduce  the  computation  load  in  subsequent  stages 
of  processing. 

Probably  the  most  representative  parameter  of  the 

phase  difference  integrity  at  a  specific  frequency  is  the 

2 

magnitude  squared  coherence  | YXy I  •  This  fact  can  be  seen 
from  Eq.  (V-35).  While  is  a  fixed  system  parameter 
|  YXy (  f n )  I  i-s  the  <^a'ta  dependent  term,  or  more  precisely 
dependent  on  the  ratio  of  the  coherent  to  uncoherent  signals 
between  the  two  sensors  considered.  Moreover  it  can  be 
related  to  the  input  SNR  under  certain  assumptions  of  equal 
and  uncorrelated  noise  between  the  sensors  and  linear  pro¬ 
cessing  up  to  the  estimation  stage,  by  the  equation  [6]: 


SNRi 


Cf  )  | 
xy  n  1 

y  (f  ) 
xy 


n 


(  V-4  6 ) 


As  such  the  magnitude  squared  coherence  ) y Xy I  can  be  a  good 
measure  for  detection  [7].  There  are  two  different  processing 
paths,  each  with  slightly  different  requirements  from  the 
detection  process: 
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a.  Time  Delay  Estimation  Path 

This  path  will  require  probably  a  lower  threshold 
on  the  coherence  because  it  is  estimating  the  time  delay  from 
a  large  number  of  frequency  bins,  as  it  will  be  shown  in  the 
next  section,  and  it  further  reduces  the  noise  level,  before 
the  trace  function  processing. 

b.  Phase  Difference  Path 

This  path  will  have  a  higher  threshold  because  it 
is  intended  to  process  mainly  line  or  narrow  band  components 
which  if  they  are  present  exhibit  a  fairly  high  coherence. 

To  preserve  the  integrity  of  the  trace  function 
which  includes  all  the  pairs  in  the  array  it  is  not  enough  to 
detect  a  potential  frequency  bin  on  a  "per  pair"  basis.  This 
operation  has  to  be  done  across  the  whole  array  (all  the  pairs). 
Two  possible  ways  to  perform  this  task  are: 

-  At  each  frequency  bin  each  pair  is  thresholded  for 
its  magnitude  squared  coherence.  Then  the  number 
of  pairs  which  pass  the  threshold  are  counted  and 
if  this  figure  exceeds  a  predetermined  number  this 
frequency  bin  is  qualified  for  further  processing. 

-  At  each  frequency  bin  the  sum  of  the  magnitude  square 
coherence  of  all  the  pairs  is  calculated  and  if  this 
sum  exceeds  a  threshold  this  frequency  bin  is  quali¬ 
fied  . 

It  should  be  noted  that  this  detection  process  based 
on  the  coherence  is  a  data  adaptive  process.  It  qualifies  the 
various  frequency  bins  based  on  its  estimated  magnitude  square 
coherence  in  real  time. 
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2 .  Compensation 

This  operation  is  straight  forward  as  explained  in 
Section  V-D. 

One  has  only  to  add  the  compensation  term  Eq.  (V-10) 
to  each  "qualified”  phase  difference  measurement  according  to 
its  frequency  and  the  participating  sensors  in  the  pair. 

3 .  Unwrapping  Problem 

This  problem  is  well  explained  in  Ref.  [32].  There 
are  several  solutions  treated  in  the  literature  [32,  36]  how¬ 
ever  in  the  case  of  the  circular  array  or  similar  two  dimen¬ 
sional  arrays  the  problem  can  be  eased  by  taking  in  consider¬ 
ation  the  following: 

a.  For  small  arrays  no  ambiguities  will  occur  at  low 

frequencies.  Ambiguities  occur  when  |t|u)>tt  for  the  largest 

aperture  (diameter  in  circular  array).  However,  for  a  small 

2  irD 

array  when  D/A<<1,  we  get  |xjw  =  — <<  it  which  is  unambiguous. 

b.  The  pairs  of  sensors  have  different  orientation 
and  one  can  find  always  (if  there  are  enough  sensors)  a  pair 
of  sensors  or  more  which  have  unambiguous  phase  difference 
measurement  even  at  higher  frequencies  than  the  limit  in  a. 
Those  will  be  the  pairs  which  are  in  a  line  parallel  to  the 
incoming  signal  wavefront. 

Using  the  above  two  facts  the  following  procedure  can 
be  outlined  to  effectively  unwrap  the  phase. 

-  Beginning  at  the  low  frequencies  where  there  is  no 
ambiguity  find  those  pair/s  which  are  perpendicular 
to  the  bearing  of  the  target/s. 
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-  At  the  higher  frequency  beginning  from  those  unam¬ 
biguous  pairs,  and  knowing  the  geometry,  i.e. 
knowing  what  is  the  biggest  possible  change  of 
phase  between  adjacent  (in  orientation)  pairs, 
one  can  readily  detect  any  2ir  jump  in  phase 
difference  and  thus  effectively  unwrap  the  phase. 

This  procedure  will  be  valid  only  if  the  target  will 
have  spectral  components  under  the  ambiguity  frequency  for 
the  given  array.  For  example,  for  a  5  meter  diameter  array 
this  ambiguity  frequency  is  150Hz  and  practically  we  can 
always  find  some  line  components  below  this  limit  (50  and 
60Hz  for  instance). 

G.  TIME  DELAY  ESTIMATION 

The  time  delay  estimation  process  is  essential  to  the 
derivation  of  the  time  delay  trace  function.  In  some  sense 
this  process  can  be  more  powerful  in  a  detection  stage  when 
the  signals  are  relatively  weak  (low  SNR  per  frequency  bin) 
because  unlike  the  phase  difference  trace  function  it 
"integrates"  over  all  "qualified"  frequencies. 

This  subject  is  extensively  treated  in  the  current 
literature  and  two  basically  different  ways  emerged: 

-  The  Generalized  Correlation  method  which  was  treated 

in  depth  by  G.  C.  Carter  [6],  several  works  derived  and 
extended  from  it  [19]  as  well  as  other  works  on  optimal 
and  suboptimal  versions  of  correlation  schemes  [8,  15]. 

-  Minimization  of  some  weighted  cost  function  [9,  14]. 


Both  the  above  methods  are  closely  related  to  the  spectral 
estimation  which  was  reviewed  in  a  previous  section  (V.E). 

During  this  work  both  methods  were  analyzed,  however,  it 
seems  that  the  second  method  is  more  applicable  to  the  problem 
of  small  arrays  mainly  because  its  time  delay  resolution  is 
not  limited  by  the  sampling  interval  as  in  the  correlation 
method.  (In  the  correlation  method  the  problem  can  also  be 
overcome  by  elaborated  interpolation  schemes . ) 

As  a  result  the  second  method  will  be  reviewed  here, 
specifically  the  development  done  by  Chan  et.  al.  [9]  which 
uses  the  weighted  least  square  error  criteria.  This  method 
was  adopted  and  will  be  followed. 

The  model  of  the  signal  to  be  processed  is  represented  by: 


x(t)  =  s(t)  +  n, (t) 

1  ( V-47 ) 

y(t)  =  as(t-x)  +  n2(t) 

The  signal  s(t)  is  uncorrelated  with  the  corrupting  noises 
n^  and  n2  which  are  stationary  random  processes. 

Assume  that  the  two  inputs  x(t)  and  y(t)  are  passed 
through  a  cross  spectral  estimation  process,  as  described 
before.  The  phase  difference  of  the  cross  spectra  and  the 
magnitude  squared  coherence  function  are  calculated.  These 
two  parameters  will  be  used  for  this  processor 


A*xy(fnJ  =  Twn+en  n  s 

V  °»FFT  •••* 

however  the  0  and  it  frequencies  will  not  be  used. 


( V-48 ) 
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en  is  an  error  (noise)  with  the  following  statistics: 

Without  signal  -  is  uniformly  distributed  from  -ir 

to  if 


2 

E  (e  }  =  0  Var{e  }  =  i- 

n  n  3 


( V-49 ) 


With  signal  for  practical 
we  can  assume  it  is  approximately 


cases  of  K  (EDF  300  ) 
gauss ian  and  has  a  variance 

(V-50) 


where 


*.  * 


n 


■wvi 


(V-51) 


The  errors  for  different  frequencies  are  uncorrelated 
i.e.  E{ej£^}=0  for  i/j  and  for  QT>8  [16], 

The  weighted  least  squares  cost  function  is: 

M 

J  =  V  (A*  -tw)2  ( V—  5  2 ) 

n  n  n 

n=l 

where  the  subscripts  of  A<}>  (f  )  were  dropped. 

xy  n 

M 

0  =  -1—  =  2  \  -Ip  (Alp  —Tit)  ) u)  ( V—  5 3 ) 

3t  n  n  n  n 

n=l 
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and 


M 

2  Wn 

n=l 

IT1 

Y]  \p  oj  2 
yn  n 

n=l 


(V-54) 


This  estimator  has  two  prominent  features: 

-  It  weighs  those  measurements  which  have  less 
variance  more  heavily. 

A 

-  It  implicitly  assumes  that  A$o  =  0  so  it  provides 
one  fixed  point  on  the  straight  line  A<f»=uit,  a 
feature  that  intuitively  should  reduce  the  variance 
of  the  estimate. 

The  statistics  of  the  estimate  are: 


/v 

E{t}  =  t  -*•  unbiased 


M 


\  '  ip  2  u>  2E{e  2} 
Tn  n  n 


E  {(  T-T  )  2  }  = 


E 

Ln=i 


^n  wm 


TT 


( V-55 ) 


K  /  <*> 

Lj  n  n 


n=l 


which  for  a  reasonable  4/  is  very  small  for  practical  values 
of  K  and  M. 

There  are  some  points  to  be  remembered: 

-  The  factor  t|i  is  not  given  but  calculated  from 


A 


4>  = 


iyi2 

i-lCl* 


( V-56 ) 
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It  can  be  defined  by  the  detection  (qualifying)  process 
mentioned  in  the  previous  section  and  as  such  the 
performance  of  the  estimator  may  be  tailored  as  de¬ 
sired. 

M 

-  The  summation  over  the  frequency  (  Z  )  does  not 

n=l 

necessarily  include  all  frequencies,  only  those  which  are 
qualified  by  the  detection  process  and  for  the  others  \p 
should  be  set  to  zero. 

This  section  will  not  be  complete  without  discussing 
the  multiple  target  problem.  For  a  multiple  target  environ¬ 
ment  the  estimator  as  defined  here  will  not  function  properly 
because  it  will  not  have  one  distinct  line  A<p  =  u>t  to  estimate. 


In  order  to  apply  the  estimator  some  preselection  of 

A 

the  A <p  measurements  should  be  done, 
n 

A  heuristic  way  to  perform  the  preselection  is 
suggested  as  follows: 

-  Beginning  from  the  lowest  "qualified"  frequency 
assume  a  "target"  line  (t^)  which  correponds  to 

A 

the  A<j>n  at  this  point  and  to  the  zero  frequency 
point. 

A 

-  Take  the  next  valid  A<j>  and  test  if 

n 


Tl"5<u,nA$n<Tl+?  ( V-57 ) 

where  C  is  chosen  based  on  some  knowledge  of  the 

a  2 
variance  of  A<J>n  i.e.  the  threshold  on  |y[  in 

the  qualifying  process. 
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-  If  it  passed  the  test  it  is  attributed  to  the  same 


"target”,  x^. 

-  If  the  test  is  not  passed  a  new  "target"  line  is 

A 

assumed  (Xj)  from  A<J>n  and  the  zero  frequency. 

A 

-  Take  the  next  valid  A<f>n  and  perform: 


tj-i  <  Un4*n<  TjH 

if  passed  attribute  to  x^ 
if  not  test  for  x^: 
t2-5  <  «.n4»n<  V? 

if  passed  attribute  to  x^ 

if  not  assume  a  new  target  line  x3 


-  When  this  procedure  is  finished  for  all  qualified 
frequency  bins  the  data  points  attributed  to  each 
"target"  line  may  be  counted  and  for  those  targets 
which  pass  a  certain  limit  of  valid  data  points 
the  estimator  in  Eq.  (V-54)  is  applied. 

This  procedure  must  be  done  for  all  the  pairs  in  the 
system  in  order  to  derive  the  time  delay  trace  function. 

It  should  be  mentioned  however  that  once  some  targets 
are  detected  the  test  will  be  performed  relative  to  them  first 
and  if  the  "new"  data  does  not  fit  within  the  limits  then  new 
"target"  lines  are  assumed. 
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H.  SYSTEM  LEVEL  PERFORMANCE 

Once  the  trace  functions  are  derived  the  bearing  estima¬ 
tors  derived  in  Chapter  III  may  be  applied.  The  resultant 
bearing  along  with  the  frequency  and  power  spectral  estimate 
for  that  frequency  can  be  displayed.  In  this  section  some 
system  level  performance  curves  will  be  presented.  System 
level  means  that  the  output  beamwidth  as  defined  in  Chapter  IV 
is  plotted  versus  the  input  SNR.  The  relaticns  used  to  reflect 
the  phase  difference  estimation  variance  back  to  the  input 
SNR^  are  Equations  (V-35)  and  (V-46)  and  their  assumptions 
set  forth: 

Eq.  (V-46) 

SNRi  =  r%- 

SNR. 

lv|  =  _ _ i_ 

|Yl  sHrT+T 

Eq.  (V-35) 

■  ttX 

SNR.2+2SNR.+1-SNR.2 

_  1 _ x _ 1 

SNR.2  Ka 
x  A 


SNRi  2 
1"(SNR.+1) 

_ l 

(SNR.  )2 

- 1 - 7  KA 

(SNRi+l)  A 


1+2SNR. 
_  x 

SNR.2  K. 
i  A 


(V-58) 
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was  chosen  at  a  value  which  will  validate  the 
gaussian  assumption  for  v(i,j)  as  defined  in  Chapter  IV,  and 
it  was  set  to  250. 

Figures  V-4  through  V-12  show  these  performance  curves 
for  circular  and  linear  arrays  along  with  the  same  CRLB  pre¬ 
sented  in  Chapter  IV  but  also  related  to  the  input  SNR. 

As  expected  these  curves  show  the  same  basic  behavior 
as  the  estimator  level  performance. 

It  should  be  remembered  that  the  relations  used, 
expecially  Eq.  (V-46)  are  somewhat  restrictive  in  their  appli¬ 
cation  and  as  a  result  the  system  performance  curves  shown 
may  not  reflect  accurately  the  real  situation.  Nevertheless 
they  can  give  a  good  perspective  on  the  expected  performance. 
NOTE:  The  SNR^  is  per  analysis  bandwidth  and  not  total  SNR 
at  the  input . 


Circular  Array-System  Level  Performance 
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Figure 


Circular  Array-System  Level  Performance 
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Figure  V-ll 


Linear  Array-System  Level  Performance 


VI.  CONCLUSION 


A.  SUMMARY  OF  THE  PRESENT  WORK 

The  aim  of  the  work  performed  in  this  dissertation  was  to 
develop  a  solution  to  the  problem  of  signal  processing, 
particularly  spatial  processing,  of  small  arrays  at  low  fre¬ 
quencies.  Special  emphasis  was  given  to  the  circular  array 
geometry  which  was  motivated  by  its  wide  practical  use. 

As  presently  available  methods  looked  unsatisfactory  to 
the  author,  a  new  concept  called  here  the  "trace  function" 
was  formulated  and  developed  in  this  work. 

It  is  based  on  the  fact  that  the  directional  information 
of  an  incoming  plane  wave  is  mainly  contained  in  the  received 
signal  phases,  or  to  be  more  precise,  in  the  relative  phase 
differences  between  received  signals  at  the  various  elements 
of  the  array. 

Those  phase  differences  are  representative  of  the  time 
delay  of  the  plane  wave  between  the  elements  of  the  array. 

The  notion  is  that  the  phase  differences/time  delays  of  the 
signals  received  from  an  array  of  sensors  which  has  a  given 
defined  geometry  describe  a  certain  "trace  function"  which  is 
constant  for  that  array  at  a  given  look  direction.  If  one  can 
find  a  geometry  for  which  the  trace  function's  basic  shape  is 
not  dependent  on  the  look  direction,  but  only  some  of  its 
parameters  are,  then  one  has  a  powerful  method  to  correlate 
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(compare)  the  trace  function  of  the  received  signal  to  a  stored 
replica. 

The  concept  of  the  trace  function  was  formalized  and  its 
application  to  some  typical  array  geometries  was  demonstrated. 
Furthermore,  it  was  shown  that  for  highly  symmetric  arrays  like 
the  circular  and  linear  arrays  the  trace  function  reduces  to  a 
particularly  simple  form. 

This  characteristic  was  used  to  derive  a  relatively  simple 
and  manageable  MMSE  estimator  for  the  bearing  of  the  in¬ 
coming  signal.  The  estimator  is  applicable  to  either  narrow 
band  signal  by  use  of  phase  difference  trace  function  or  to 
the  wide  band  signal  using  time  delay  trace  function. 

The  performance  of  the  estimator  was  checked  by  simulation 
and  compared  to  the  CRLB  as  adapted  to  this  application. 

A  brief  performance  summary  is: 

-  The  estimator  exhibits  a  threshold  phenomena  called 
"lock-on"  with  respect  to  the  variance  of  the  phase 
difference  measurements. 

-  Within  the  "lock-on"  range  the  estimator  is  practically 
unbiased. 

-  With  respect  to  beamwidth  performance  the  estimator  is 
asymmtotically  efficient. 

-  It  exhibits  a  saturation  point  in  the  number  of  elements 
at  constant  diameter  or  in  the  diameter  at  constant 
number  of  element,  beyond  which  only  marginal  improve¬ 
ment  is  achieved  in  beamwidth. 


The  estimator  was  applied  to  both  circular  and  linear 
arrays  and  its  performance  was  equally  good  except  for 
differences  imposed  by  the  geometry.  (symmetry) 

Finally,  a  system  configuration  applying  trace  function 
principles  was  outlined  and  the  major  problem  areas  caused  by 
the  specific  application  were  identified,  reviewed  and  some 
suggestions  to  solutions  were  made. 

As  a  final  conclusion  it  can  be  said  that  this  work  pre¬ 
sented  a  new  way  of  looking  at  the  available  information  from 
an  array  of  sensors.  It  is  clear  that  the  use  of  the  above 
concept  in  a  generalized  replica  correlation  model  is  effective 
for  highly  symmetrical  arrays.  As  such  the  circular  array  is 
a  natural  one  to  be  treated  this  way.  For  this  array  the 
derivations  of  the  bearing  estimator  is  simple  and  its  per¬ 
formance  appears  to  be  very  good. 

However,  it  should  be  emphasized  that  the  trace  function 
concept  is  applicable  to  any  array  geometry  and  for  both 
narrow  and  wideband  processing  schemes. 

B.  OTHER  POSSIBLE  APPLICATIONS 

The  trace  function  concept  and  the  estimators  derived  in 
this  work  are  not  limited  in  their  application  to  the  passive 
sonar  case. 

They  can  be  applied  to  related  areas  like  seismic,  and 
ocean  wave  analysis  as  well. 


In  addition  the  principle  can  be  applied  to  any  other 
situation  in  which  phase  differences  and/or  time  delays 
within  an  array  of  sensor  can  be  measured.  The  performance 
curves  show  that  the  trace  function  bearing  estimator  is 
very  tolerant  to  phase  difference  measurement  noise,  even 
30°-40°  of  standard  deviation  in  the  measurements  can  give 
very  good  results  with  reasonable  number  of  elements.  This 
suggests  that  even  in  areas  where  F.F.T.  techniques  for  cross 
spectral  estimation  are  still  beyond  the  state-of-the-art 
(like  microwave  frequencies)  available  hardware  exists  for 
phase  difference  or  time  of  arrival  measurements  which  may 
be  sufficient  to  allow  application  of  the  trace  function 
concept . 

C.  PROPOSED  FUTURE  WORK 

Although  the  trace  function  principle  was  developed  in 
this  work  there  are  a  vast  number  of  points  to  be  checked. 

In  addition  several  extensions  and  applications  should  be 
investigated . 

Some  proposed  topics  are: 

1.  Testing  the  concept  of  3  possible  ways: 

a.  Computer  simulation  of  a  whole  system 
(similar  to  Figure  V-l)  including  model  of 
the  sea. 

b.  Laboratory  simulation. 

c.  Testing  on  real  data  at  sea. 

2.  Investigation  of  the  weighted  MMSE  ,  optimi¬ 
zation  of  the  weights  and  their  possible  adaptive 
derivation. 

. _  _  _ _  .  'i. _  . 
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3 .  Derivation  of  a  minimization  algorithm  for  "non- 
analytic"  trace  function. 

4.  Further  formalization  and  development  of  the 
phase  preserving  spatial  decorrelation. 

5.  Investigation  of  the  multiple  target  discrimina¬ 
tion  for  both  narrow  and  wide  band  targets . 

The  above  are  only  a  few  of  the  possible  topics  opened  up 
by  the  introduction  of  the  trace  function  principle. 


APPENDIX  A 


TRACE  FUNCTION  DERIVATIONS 
1.  Circular  Array 

This  derivation  is  related  to  Figure  II-4.  Using  Eq.  CII-2) 
the  following  is  the  interpretation  of  the  parameters  for  cir¬ 
cular  arrays : 

A  A 

ci  »  -sinQx  +  cos0y 

r^  =  R  C-sin9^p  ix  +  cos  ^p  iy] 


A  A 

where  x,  y  are  unit  vectors  in  x  and  y  directions. 
Inserting  in  Equation  (II-6) t 


At..  =  p  (-sin8x 

1 J  U 


+  cos9y)  •  [-(sin-|p  I  -  sin-^p  j)  x 


,  2  TT  • 
(cos-jf  1 


2  IT  .  . 
-  COS-JJ-  3) 


y  ] 


a  ^  Ccos0(cos^p  i  -  cos-^p  j)  +  sin9(sin-pp  i  -  sin-^p  3] 
=  |  CcosC^p  I  -  9)  -  cost—  j  -  9)] 

=  -2§  sinCj  (i+j)  -  9]  sin[|  (i-j)]  (A-l) 
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and  for  phase  difference  Eq.  (II-5): 

*  4Tij  “i  =  -Tf  sin  CJf  Ci+j)  '  e]sinCf  < 

2 .  Linear  Array 

This  derivation  is  related  to  Figure  11-15.  The 
are  the  various  terms : 

A  A 

a  =  sin0x  +  cos8y 

A 

r .  =  d  i  x 

i 

Inserting  in  Eq.  (II-6): 

d  A  A  A 

At..  =  —  (sin0x  +  cos0y)*[(i  -  j)x] 
i]  c 

=  -  sin0(i-j ) 
c 


and  for  phase  difference  Eq.  (II-5): 


A$..  =  At.,  u  =  3^  (i-j)sin0 


-j ) 3  (A- 2 


following 


( A-  3 ) 


( A-4 ) 
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APPENDIX  B 


1 .  Computer  Program  for  Three  Dimensional  Plotting  of  the 
Trace  Functions  for  Circular  Arrays. 

The  following  computer  program  as  well  as  all  the  other 
programs  in  this  dissertation  were  developed  and  executed  on 
a  Hewlett-Packard  9845T  computer. 
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THREE  DIHENSIOMBL  TRfiCE  FUNCTION  PLOTTING  FOR  CIRCULftR  fiRRNV 
PARAMETER  DEFINITION 
H  -  Number  Of  Elements 
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2 .  Computer  Program  for  Three  Dimensional  Plotting  of  the 
Trace  Function  for  Linear  Arrays. 

Following  is  the  listing  of  the  computer  program  for  the 
plotting  of  the  trace  function  for  linear  arrays. 
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3.  Computer  Program  for  Three  Dimensional  Plotting  of  the 
Trace  Function  for  Random  Arrays . 

Following  is  the  listing  of  the  computer  program  for  the 
plotting  of  the  trace  function  for  random  arrays.  The  same 
program  may  be  used  for  conformal  arrays . 
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APPENDIX  C 


SIMULATION  AND  THE  NOISE  MODEL 


The  simulation  of  the  performance  was  done  by  generating 
a  noisy  trace  function  (Eq.  III-l). 


A<f>ij  =  fTr(i,  j  ,6)  +  v(i,  j  ) 


(C-l) 


and  applying  the  MMSE  estimator. 

The  trace  function  fTr(i,j,0)  is  generated  simply  by 
calculating  its  value  for  all  " i"  and  "j"  at  a  given  bearing 
”0".  The  additive  noise  v(i,j)  is  the  more  problematic  term 
and  its  definitions  will  be  discussed  next. 


1.  Noise  Model 

In  Ref.  [29]  it  was  shown  that  for  high  SNR  or  conversely 
for  low  SNR  with  long  averaging  time  (high  K)  the  error  in  the 
phase  difference  estimate  may  be  considered  approximately  as 
a  gaussian  distributed  noise  with  zero  mean  and  the  variance 
given  by  Eq.  (V-35) 


>iy‘ 

ka 


(  C-2) 


at  a  specific  frequency. 

This  assumption  of  gaussian  noise  was  used  throughout  the 
simulation  and  to  validate  it  was  chosen  to  be  250. 
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1 


This  K  corresponds,  by  using  Eq.  (V-37),  to  a  T  (data 
stationarity  time)  of  approximately  85  seconds  when  assuming 
one  hertz  frequency  resolution  in  processing.  This  figure 
is  a  reasonable  one  for  real  data  in  the  ocean. 

_3 

The  values  of  coherence  square  used  were  in  range  ~  10 
to  0.985  which  corresponds  to  input  SNR  of  -15  to  +21  dB 
(using  Eq.  (V-46)). 

The  next  point  to  be  considered  is  the  covariance  between 
the  various  noise  terms .  In  the  phase  difference  estimation 
process  each  element  is  used  in  2N-1  pairs  so  that  obviously 
there  is  a  certain  correlation  between  these  pairs.  The 
following  model  is  assumed  for  the  covariance  matrix  which 
takes  in  consideration  the  various  pairs  of  terms  of  the 
trace  function.  Table  C-l  summarizes  the  above  model. 


r— 

! 

DESCRIPTION 

TERM  FORM 

COVARIANCE 

a 

Pairs  without  common  elements 

A  A. 

Ad> .  .  and  A<f>.  , 

Yi]  kH 

0 

2 

Pairs  with  one  common  element 

A$  .  .  and  A$  .. 

Tik 

A$^j  and  A$k^ 

0.5 

-0.5 

3 

Pairs  with  two  common  elements 

Acb .  .  and  A$  .  . 

A$ . .  and  A$ . . 

ID  33 

0 

0 

1 

Pairs  with  all  common  elements 

A  A 

Ad> .  .  and  Ad>  •  • 

13  31 

A$£^  and  A$^ 

-1 

1 

Table  C-l 


The  covariance  array  is  basically  a  four  dimensional  array 


which  takes  in  consideration  all  the  above  pairs  of  terms. 

But  for  ease  of  perception  and  calculation  this  array  was 
described  as  a  two  dimensional  array  which  is  composed  of 
submatrices . 

.  2  2 
This  covariance  matrix  is  of  order  N  xN  and  corresponds 

2 

to  the  N  noise  terms  of  the  trace  function. 

An  example  of  such  a  covariance  matrix  for  a  4  element 

array  (16  term  trace  function)  is  shown  in  Figure  C-l. 

The  method  in  Ref.  [4]  was  used  to  obtain  correlated 

gaussian  random  variables  having  a  specified  covariance  matrix. 

2  2 

Specifically,  defining  the  above  N  xN  covariance  matrix  as  Q 
then 

Q  =  U  A  U*  (C-3) 

where  U  -  is  a  matrix  whos  •>.  columns  are  the  eigenvectors  of  Q 
A  -  is  a  diagonal  matrix  whose  elements  on  the  diagonal 
are  the  eighenvalues  of  Q 
U^-  transpose 

2 

The  next  step  will  be  to  generate  N  independent  gaussian 
random  variables  with  zero  mean  and  unity  variance . 

This  operation  can  be  done  by  several  methods,  in  the 
present  case  the  polar  transformation  method  was  used  [4]  to 
generate  a  pair  of  gaussian  independent  random  variables  from 
a  pair  of  uniform  random  variables  which  are  available  as  a 
built-in  function  in  the  computer. 
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Trace  Function  Noise  Covariance  Matrix 


To  generate  the  dependent  gaussian  variables: 

V  =  U  A5  Z  CC-4) 

where  V  -  dependent  random  variable  vector 

Z  -  independent  random  variable  vector 
As  a  result  the  covariance  of  V  is : 

E  {vv1}  =  v  a1  fitzz^A*  (c-5 ) 

E{ZZt}  =  I  -  independent,  unity  variance 

E{VVt}  =  U  A  Ut  =  Q  CC-6) 

To  get  the  required  noise  power  V  will  be  multiplied  by  the 
scalar  vffip  (noise  r.m.s.  amplitude): 

No=  vJTp  U  it1  Z  (C-7) 

2 .  Simulation 

Once  the  additive  noise  was  generated  it  was  added  to  the 
trace  function  calculated  for  the  circular  or  linear  array 
Eq.  (II-8)  and  (11-10)  respectively. 

Next  the  corresponding  MMSE  estimator  was  applied  to 
Eq.  (III-23)  and  CIII-24). 

The  above  operations  were  run  48  times  and  the  mean  and 
variance  of  the  resultant  bearing  estimates  were  calculated 
and  printed  out.  In  addition  a  plot  of  the  beamwidth  versus 
input  SNR  was  done. 
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The  CRLB  Eq.  (IV-3)  and  (IV-4)  respectively  are  given 
only  with  the  simulation  results  for  comparison  purposes. 

Following  is  a  listing  of  the  simulation  program  for  an 
8  element  circular  array.  In  order  to  change  the  number  of 
elements  in  addition  to  change  of  the  parameter  N  the  rele¬ 
vant  arrays  in  the  main  program  have  to  be  redimensioned, 
(line  240) 

In  order  to  change  to  a  linear  array,  equations  have  to 
be  changed  to  the  applicable  ones  in  the  generation  of  the 
signals,  the  estimators  and  the  CRLB  (lines  200,  360,  680, 
1270,  1670,  1680,  1710,  1760,  1770,  1800). 

The  transformation  matrix  for  the  generation  of  the 
correlated  random  variable  (t  =  U  A**)  has  to  be  calculated 
only  once  for  a  given  number  of  elements,  then  this  part  of 
the  program  can  be  skipped  (lines  260  to  290).  Care  should 
be  taken  to  have  a  data  file  (TRAN)  defined  for  the  storage 
of  the  transformation  matrix.  The  subroutine  "Symqr"  (line 
2170)  for  the  eigen-analysis  can  be  found  in  the  Hewlett- 
Packard  Numerical  Analysis  Software  Package. 
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PHASE  DIFFERENCE  TRACE  FUNCTION  ESTIMATION  -  SIMULATION  PROGRAM 
PARAMETER  DEFINITION 
D  -  Array  Diameter 
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2140  NEXT 


2150  PRINT  “*******  COVARIANCE  MATRIX  **** 
2160  MAT  PRINT  C 

21 70  CALL  Sy»qr<C<*>f D< * > , E<*> , 0, N' 2, IE-20 
2180  PRINT  "************  EIGENVALUES  ***** 
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